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SUMMARY 

This paper presents a computer program for maximum likelihood and 
Bayesian estimation of the vector p of multinomial cell probabilities 
from incomplete data . Also included is coding to calculate exact and 
approximate elements of the posterior mean and covariance matrices. 

The program is written in FORTRAN IV language for the Control Data 
CYBER 170 series digital computer system with network operating system 
(NOS) 1.1. The program requires approximately 44000 octal locations 
of core storage. A typical case requires from 72 seconds to 92 seconds 
on CYBER 175 depending on the value of the prior parameter. 
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DESCRIPTION 


This paper describes the main computer program used in reference 
1 for estimating the vector p of multinomial cell probabilities from 
incomplete data. The data is incomplete in that it contains partially 
classified observations. Each such partially classified observation is 
observed to fall in one of two or more selected categories but is not 
classified further. The estimation criterion is minimization of risk 
for quadratic loss L(p-p) = {p-p)'(p-p) for p an estimator of p. 

In addition, elements of the posterior mean and covariance matrices 
are calculated exactly and approximately. A Taylor-series function is 
used to approximate the posterior covariance matrix. A Taylor-series 
function, the maximum-likelihood estimate, and the posterior mode are 
used to approximate the posterior mean. 

Monte-Carlo simulation studies are performed for small- and medium- 
size samples to assess 

(1) which of the maximum-likelihood estimate, posterior mode, and 
Taylor-series approximate posterior mean best minimizes risk 
for specified values of p; 

(2) how well each of these functions approximates the exact pos- 
terior mean; and 

(3) how well a Taylor-series function approximates elements of 
the posterior covariance matrix. 

Samples are of size 25 and 50, percentage of incomplete data varies 
around 15 and 40, and probabilities range from the center of the prob- 
ability simplex ?2 to one of its corners. Probabilities equal the means 
of the prior distributions for varying parameters or are randomly gen- 
erated from these distributions. An exploratory robustness study is 
conducted by using the correct prior, a uniform prior, and a perturbed 
prior in the Bayesian estimators. The iterative algorithm of Dempster, 
Laird, and Rubin (ref. 2) is used to evaluate all three estimators. 
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Other discussion, analysis, and results are given in reference 1. 
Included in the discussion in reference 1 are descriptions of 
pseudorandom-number generators for the Dirichlet, uniform, and trinomial 
distributions. Also given are tree diagrams (figs. 5.1 - 5.3 of ref. 1) 
that illustrate the flow of the computer program. 

The computer used is a Control Data Corporation (CDC) CYBER 170 
series digital computer system with network operating system (NOS) 1.1. 
This computer operates with a 60-bit word and single-precision accuracy 
of about 14.5 significant figures. The programing language is FORTRAN 
Extended, Version 4.6. The program requires approximately 44000 octal 
locations of core storage. A typical case requires from 72 seconds to 
92 seconds depending on the prior parameter. 


A listing of the main computer program is given in the third sec- 
tion. An index precedes the listing. Symbols are defined in the list- 
ing. Note that the program is written for DESIGN 2 of reference 1 but 
can be modified to be DESIGN 1 of reference 1 by deleting the Dirichlet- 
generation level and changing the dimensioning of the QLMSIJ matrix. 

DESIGN 1 has a fixed-effects model constituting full factorials (4x2 x3) 

★ 

having four levels and two replications per cell. DESIGN 2 has a 
mixed-effects model constituting nested factorials (4-10x2 x3) having 
four levels and two replications per cell within each of four variations 
-of- the prior parameter y._Th^e ten generations of the Dirichlet prob- 
ability p in DESIGN 2 are considered random; the remaining factors a>¥ ' 
considered fixed. 


★ 

level 1 
level 2 
level 3 
level 4 


four Dirichlet probabilities p 
two sample sizes 

two percentages of incomplete data 
three estimators 
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Tq run a case, any necessary changes are made to the following 
lines in MAIN 

XNUl = 

XNU2 = 

IP = 

SEED = 

In addition, the subroutine GAM, which is a function of the prior 
parameter v (NU), and the Hollerith labels 

DATA ALABEL/IOHA. NU = ( ,10H0. 1,0. 1 ,9. ,10H8) . / 

DATA TLABEL/IOHTABLE 7.1 / 

are changed as needed. For values of y that are less than 10, GAM(y+l) 
is calculated from the relationship GAM(y+l) = y GAM(y) and a starting 
value. For y an integer, the starting value isGAM(2) = 1 and for 
y = X 1/3, for x an integer, the starting value is GAM(3 1/3) = 
2.7781584804296. For y greater than or equal to 10, Stirling's formula 
is used to approximate the gamma function to 11 significant figures of 
accuracy. 

Because a MODIFY system (ref. 3) is used to maintain the program 
on a permanent file, a new case is easier to make by changing lines of 
code rather than reading data cards. Outputs consist of printouts and 
tapes. Some tapes are directly used as tables. Tapes are also usually 
input to canned programs for calculating analyses of variance and to a 
program for surming biases or mean squared error over replication, sam- 
ple size, percentage of incomplete data, and/or generated Dirichlet 
probability. 

Subroutines URAN, URANV, and MATINV shown in the coding are from 
the NASA, Langley Research Center, mathematics computer library. They 
are described in Appendices A, B, and C, respectively. Subroutine URAN 
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gives a single uniform random number according to the algorithm describ-^ 
ed in Appendix A. Subroutine URANV gives a vector of uniform random 
numbers from URAN. Subroutine MATINV solves a system of simultaneous 
linear equations. 

In addition, other computer programs not given in the listing have 
been written. Among these are programs to test the gamma, Dirichlet, 
trinomial, and uniform pseudorandom-number generators; to calculate 
analyses of variance; and to sum mean squared errors and biases over 
one or more of replication, percentage of incomplete data, sample size, 
and generated Dirichlet probability. The sums have been used for plots 
in reference 1. Note that a number of subroutines from IMSL (Interna- 
tional Mathematical and Statistical Libraries, Inc.; ref. 4) have been 
used in calculating the analyses of variance. 
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INDEX OF COMPUTER PROGRAM 


NAME 

MAIN 

GAMMA 

GENXZ 

EPM 

GAM 


METHODS 


COUNTS 

ESTMSE 

KTITER 

SUMMARY 

BESTEST 


USAGE PAGE 

main program 6 

subroutine to generate a gamma random variable 35 

subroutine to generate trinomial complete (x) and 37 

incomplete (z) data; also calculates complete-data 
maxi mum-likeli hood estimate 

subroutine to calculate exact posterior mean and 41 

covariance matrices 

subroutine to evaluate the Gamma function by exact 45 

values or Stirling's approximation; note that this 
subroutine differs for each of the four different 
values of the prior parameter \>eNU (only the sub- 
routine for v=(0. 1,0. 1,9.8) is shown) 

subroutine to calculate the Taylor-series approxima- 47 


tions APM and APC for the posterior mean and covar- 
iance matrices; also calculates the posterior mode 
PMD and the maximum-likelihood estimate MLE from in- 
complete data 

subroutine for covariance approximations and complete- 54 
data maxi mum- likelihood estimate to count the number 
of the 200 trinomial simulation trials that have nega- 
tive, zero, and positive error and that have absolute 
relative difference less than certain percentages 

subroutine to calculate the usual, control -variate, 56 

and regression estimates of mean squared error and 
their sample variances 

subroutine to increment counters for averaging the 58 

number of iterations an estimator requires and for 
determining how many of the 200 trinomial simulations 
an estimator requires a specified number of iterations 

subroutine to calculate Tukey's five-point data sum- 60 
mary: median, hinges, and extremes 

by two different criteria (summed absolute relative 61 

difference and sunmed squared error), subroutine de- 
termines which estimator is best for a given one of 
the 200 trinomial simulations; also includes a section 
corresponding to COUNTS for the estimators APM, PMD, 
and MLE 
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LISTING OF COMPUTER PROGRAM 


PgOGRAN HAIN( INPUT* aUTPUT*TAPE5»TAPE12) 


PR06RANED BY KAREN RACKLEY CREDEUR* SPRING 1977* CDC 6600* 

FORTRAN EXTENDED VERSION 4.6* NASA* LANGLEY RESEARCH CENTER 

INTEGER COVSKIP 
REAL NSE(6*7) 

DIHENSIQN BESTQL<4»2*3)*CTSEQL(3) 

DINENSION R(2»7*2)*rUKEY(10)*UUt50) 

OINENSION 0LNS1(2»2*2)>QLNS2(2*2*2)»QLNS3(2*2*2)»QLHS4(2*2*2) 

OlNENSION QLNS5(2*2*21*QLHS6(2*2*2l»E6IASlt2*2*2l*EBIAS2C2*2*2) 

OINENSION EBIAS3(2»2>2)*EnSll2*2*2>*ENS2(2*2*2)*EHS3(2*2*2l 
OINENSION QLHSll(10»2*2)*QLNS12fl0*2*2)*QLHS21(10»2*2) 

OINENSION OLMS31(10*2*2)*QLNS32llO*2*2)*OLNS22llO*2*2) 

OINENSION 0LNS4K 10*2«2)*QLHS42(10*2>2I*0LHS51(10»2*2> 

OINENSION QLNS6m0»2*2l*QLNS62ll0*2*2)*OLNS52(10*2*2) 

OINENSION T11(2*2*1U>T12(2*2*11)»T21(2*2*11)*T22(2*2*11) 

OINENSION T31(2*2»ll)*T32<2*2*ll)*T41(2»2»in*T42(2*2»ll) 

OINENSION T5l(2*2»ll)»T52(2*2*ll)*T61(2*2*lll*T62(2»2*ll) 

OINENSION ALABEL(3)*TLABEL(1) 

CONNON OEP(3*3)»OOL(4»3)*E(4*3I»IROBUST*NTS*P1»P2*P3*TPID 
C0NN0N/BEST/BESTEP(3»2)»CTR0EP*CTR0QLC3)*PRDEP<9*3)*PR0QL(9*7I*SBI 
1ASEP(3*3I*SBIASQL(3*7) 

CONNON/BlASRO/COUNTB(3»e)»COUNTRO(8*8) 

C0NN0N/CALEST/APNC11*APNC12* APNC22*C0NVCRI*C0VSKIP*0NLC1*DNLC2»DNL 
lC3*0PID*EAPMCll*EAPMC12*EAPNC22»EPNCll*EPMCl2*EPNC22*IST0P»Ni2»N13 
2»N23»PI0*PNLCl»PMLC2»PMLC3*SS*SSN»TINAP»TIHEP*TIH<2»*TIH2l*TIH31*X 
3NU1*XNU2*XNU3»Z1*Z2*Z3*Z12*Z13*Z23*Z1N*Z2N*Z3N 
C0NN0N/ITKT/AVNUNIT(6)»CTNUNIT(6»10) 

EOUI VALENCE ( E ( 1* 1) * PEPNl I* (E (1*2 ) * PEPN2 )* ( E (1* 3 )» PEPN3 ) 

EQUIVALENCE ( E ( 2* 1 ) * PNLE 1 1* ( E ( 2* 2 1 * PNLE2 ) * ( E (2* 3 ) * PNL E3 ) 

EQUIVALENCE ( E ( 3* 1 ) » PPNDl )* ( E (3* 2 ) * PPN02 >* ( E ( 3 * 3 > * PPN03 ) 

EQUIVALENCE (E(4*l)iPAPNllV(E(4y2)»PAPN2»* (E (4*31* PAPH3) 

EQUIVALENCE ( DEP (1> 1 ) » ENLl) * (DEP(1*2)*ENL2)* ( DEP (1 * 3 ) * ENL3 ) 

EQUIVALENCE ( DEP ( 2* 1 )* EPNOl )* ( DEP ( 2* 2 ) * EPNDZ ) » (DEP (2 * 3 ) * EPN03 ) 

EQUIVALENCE (DEP ( 3* 1 ) * E APNNl )* ( DEP ( 3*2 ) * E APNN2 > » ( OEP( 3* 3 )» EAPNN3 ) 

EQUIVALENCE ( OQL ( 1* 1 )* DE PNNI }*( DOL ( 1*2 )* DEPNN2 >» ( DQL ( 1* 3 1 * DEPNN3) 

EQUIVALENCE ( DQL ( 2* 1 ) * ONLl) * ( DQL ( 2 *2 ) * DNL2 ) * ( DQL ( 2* 3 ) *DNL 3) 

EQUIVALENCE ( DQL ( 3* 1 ) * OPNDl) * ( DDL ( 3*2 ) * DPND2 1 * (DQL ( 3* 3 ) * 0PH03 ) 

EQUIVALENCE (DQL ( 4* 1 )* DAPNNl I* (DQL (4*2 ) *DAPNN2 ) * ( DQL ( 4* 3 ) »0APNN3 t 

DATA ALABEL/IOHA. NU • (* lOHO. 1* 0 .1* 9. * 10H8) / 

DATA TLABEL/IOHTABLE 7.1 / 
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CONVCRl RELATIVE'ERROft CONVERGENCE CRITERION (USUALLY 0.0001) 

IP DENOTES. IN THE FOLLOWING ORDER. ONE OF THE EXPECTED P"S 

(.01..01..98)»(.1..1..6).<.2>.3*.5). AND ( 1 /3. 1 /3. 1/ 3 ) 
NSS INTEGER SS 

NXZ NUMBER OF TRINOMIAL (X AND Z DATA) SIMULATIONS 

PAPMI I-TH T.S. APPROXIMATED POSTERIOR MEAN 

PEPHI I-TH EXACT POSTERIOR MEAN 

PI I-TH GENERATED P 

PIO PERCENTAGE OF INCOMPLETE DATA 

PMLCI I-TH COMPLETE M.L.E. 

PMLEI I-TH INCOMPLETE-DATA M.L.E. 

PPMDI I-TH POSTERIOR MODE 

SS SAMPLE SIZE 

SSN SS ♦ SUM OF PRIOR PARAMETERS XNUI 

XNU VECTOR OF PRIOR PARAMETERS XNU> ( XNUli XNU2. XNU3 ) 

ZI NUMBER OF OBSERVATIONS FALLING IN CATEGORY I 

ZIJ NUMBER OF OBSERVATIONS SUCH THAT EACH OBSERVATION IS 

KNOWN TO FALL IN ONE OF CATEGORIES I AND J BUT IS NOT 
FURTHER CLASSIFIED 
ZIN ZI«XNU1 


XNU1*0.1 

XNU2-0.1 , ‘ 

IP-2 

SEED-24158739. i 

GSEED-SEED+100. . 

INITIALIZE ONE-OIMENSIONAL FORM OF UNIFORM RANDOM-NUMBER 
GENERATOR FOR GENERATING GAMMA RANDOM VARIABLES 

UN-URAN(GSEED) 

PRINT 4i GSEED.UN 

XNU3-10.-XNU1-XNU2 

XNU-XNUl^XNUZi-XNUB 

GENERATE A 3-COMPONENT (2-DIM) VECTOR OF DIRICHLET PROBABILITIES 

DO 9910 IGEN-1.10 
2 Gl-GAMMA(XNUl) 

G2-GAMMA(XNU2) 

G3-GAMMA(XNU3) 

G-61«G24>G3 
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P1>G1/G 
P2-G2/G 
P3-1.-P1-P2 
IF (P3-1.) 7»7.5 

5 PRINT 3« XNU1»XNU2*XNU3>XNU«IP«P1»P2»P3>PID>NSS*NXZ»C0NVCRI 
PRINT 6 

6 FORMAT!* P3 IS NEGATIVE. REGENERATE DIRICHLET.*//// ) 

GO TO 2 

7 PIO-0.15 
NSS-25 
NXZ-200 

CONVCRI-0.0001 
KASE-0 
IPIO-1 
ISS-l 
IPRINT-0 
1 SEEO-SEED*2« 

PRINT 3* XNU1»XNU2>XNU3»XNU*IP»P1*P2>P3»PID»NSS»NXZ»C0NVCRI 
3 FORHATUHl** XNU1-*F6.3* XNU2«*F6.3* XNU3«*F6.3* XNU»*F6.2* IP»*I2 
I* P1-PF6.A* P2»*F6.^* P3-PF6.A* PI0»*FA.2* NSS»*I3* NXZ-*I3* CONVC 
2RI»*F7.5//I 

INITIALIZE VECTOR FORM OF UNIFORM RANDOM-NUMBER GENERATOR 

CALL URANV(SEEO#l»UN) 

PRINT SEEOfUN 

A FORMAT!* SEED GIVEN URANV IS*E23.1A* SEED TRANSFORMED BY URAN FRO 
IM THIS SEED IS*E23.1A7/) 

SS-NSS 

SSN«SS*XNU 

DO 95 NREPLIC«1»2 


NUMXZ«NXZRl«NXZR2-NCOV-NXZ 

ISTOP-0 

COVSKIP-0 

AVTPID«AVOPID«0. 

INITIALIZE COUNTERS FOR X RELATIVE DIFFERENCE AND SIGN OF BIAS 

DO 28 I-l»9 
DO 27 J-l»7 
PRDQLII* J)*0. 

27 CONTINUE 

SBIASEPm-0. 
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28 CONTINUE 

DO 29 1-1*21 
SBlAS0L(I)-0. 

29 CONTINUE 

00 30 1-1*27 
PR0EP(I)-0. 

30 CONTINUE 

CTR0EP-CTft0QL(li-CTR0OL(2t-CTR0QL(3)-0. 


INITIALIZE COUNTERS FOR NUNBER OF ITERATIONS FOR CONVERGENCE 

00 ^0 11-1*6 
AVNUHIT(Il)-0. 

00 39 12-1*7 

39 CTNUHIT(Il«I2l-0. 

CTNUNITm*8)-CTNUMIT(Il*9)-CTNUMlTni*10»-0, 

40 CONTINUE 

TINm-TIM(2>-TlNEP-TIHAP-0. 

00 42 J-l*2 

BESTEP(l*J)-BESrQLIl* Jf l)-BE$T0L(l»J»2)*BEST0L(l»J*3)-0. 
BESTEP<2* J}-BESTQL(2*J*1)-BEST0L(2* J*2)-BESTQL(2*J*3)-0. 
BESTEP(3*J)-BESTOLC3*J*1)-BGSTOLO»J»2)-6ESTQL(3*J*3)-0» 
BESTQL(4*J*1)-BESTQL(4*J*2)-BEST0LI4* J*31-0. 

42 CONTINUE 

INITIALIZE COUNTERS FOR BIAS AND RELATIV.E DIFFERENCES 

DO 49 K2-l»8 
00 44 Kl-1*8 
C0UNTR0<Kl*K2)-0. 

IF CKl-3) 43*43*44 

43 C0UNTB(Kl*K2l-0. 

44 CONTINUE 

45 CONTINUE 

INITIALIZE AVERAGES AND ERROR SUMMARIES TO ZERO 

AVERAGE ESTIMATORS AND THEIR STANDARD ERRORS fS.E.) 

APHLEl- APEPNl- APAPMl- APPMDl- APMLCl- APl-0. 
XAPHLEl-XAPEPMl-XAPAPMl-XAPPMOl-XAPMLCl-XAPl-O. 

APMLE2- APEPM2- APAPM2- APPM02- APMLC2- AP2-0. 
XAPNLE2-XAPEPN2-XAPAPM2-XAPPMD2-XAPHLC2-XAP2-0. 


s. 
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AVPMN1*AVPNN2«AQPHN1«AOPHN2>AOI>H01>AOPHD2*0. 
AFPNN2«AFPNN2«A6PMN1*AGPHN2> PGN01« PGN02«0. 

AV BIAS* nSE» e VAR(HSE) OF ESTIMATORS FROM EXACT POSTERIOR MEAN 

XAPNN1«XAPHN2>XAPMNA>XAPHNSQ<0. 

XPNOl «XPM02 «XPMOA -XPHOSO «0. 

XNLl *XHL2 aXMLA -XMLSQ -O. 

AV BIAS* MSE* C VARfHSEl OF ESTIMATORS FROM GIVEN OR GENERATED P 

YEPMNl -YAPMNl «YPM01 -YMLl -YMCCl «VAPMNl -OARMNl -OPMOl «0. 

YEPMN2 •YAPMN2 -YPM02 -YML2 *YNLC2 •VAPMN2 •OAPMN2 ■OPHD2 -O. 

YEPMNA -YAPMNA -YPHOA ■YMLA -VAPHNA "OAPHNA -OPMOA -0. 

YEPMSQ -YAPMSO -YPHOSO-YMtSO -VAPMSQ -OAPMSO ■OPMOSO-0. 

RE6APMO-REGPMOO*REGEPMO«REGMLO> REGMLC0*YHLC20«0. 

REGAPM1«R£GPM02> REGMLCl" YHLC21-0. 

REGAPM2- REGMLC2*YMLC22*0, 

INITIALIZE FOR S.E. OF AVERAGE BIAS 

WAPMN1«WPM01>WML1> WAPMN2-UPM02-WNL2-0. 

UNLCI»UEPHN1«UAPHN1»UPM01-UMLI-FAPMNI-GAPMM1«GPM01»0. 
UMLC2«UEPMN2-UAPMN2*UPM02>UML2«FAPMN2-GAPMN2«GPMD2*0. 

EST AVERAGE* BIAS* S.E.* AV t REL DIFF* ( MSE FOR COV EST FROM EP 

AEPHCll-AAPMCll»BAPMCll-CllMSE-VCllMSE-APROCll-0. 

AEPMC12-AAPMC12-BAPMC12-C12MSE-VC12MSE-APRDC12-0, 

AEPMC22-AAPMC22-BAPMC22-C22MSE-VC22MSE-APRDC22-0. 

A2EPCII-A2EPC12-A2EPC22-VARMSE-0. 

A2APCI1-A2APC12-A2APC22- VARA-0. 


GENERATE NUMXZ SETS OF TRINOMIAL DATA* CALCULATE ESTIMATORS* AND 
COMPARE THEM AS APPROXIMATIONS FOR THE EXACT POSTERIOR MEAN. 

ALSO ASSESS HOW WELL THE T.S. APPROXIMATION DOES FOR THE EXACT 
POSTERIOR COVARIANCE MATRIX. COMPARE ESTIMATORS WITH GENERATED 
P AND DO EXPLORATORY ROBUSTNESS STUDY 


DO 65 NT>1*NXZ 
NTS-NT 
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GENERATE NULTINOMIAL COMPLETE DATA X AND INCOMPLETE DATA Z 
CALCULATE COMPLETE>OATA MAXIMUM LIKELIHOOD ESTIMATE 

CALL 6ENXZ(UU*NSS) 

If ClSTOP-11 46*58,58 

CALCULATE EXACT POSTERIOR MEAN 

46 CALL EPM 

fOR INCOMPLETE DATA Z* CALCULATE MAX IMUM>LIKEL IHOOD ESTIMATE# 
POSTERIOR MODE# AND TAYLOR-SERI ES APPROXIMATED POSTERIOR MEAN 
AND TAYLOR-SERIES APPROXIMATED POSTERIOR COVARIANCE MATRICES 

CALL METHODS 

IFIISTOP-II 49,58,58 

AVERAGE ESTIMATORS 

49 APMLEI»APMLE1*PMLE1 
APMLEZ«APMLEZ+PMLE2 
APEPHl-APEPHl^PEPMl 
APEPM2«APEPM2*PEPMZ 
APAPMl-APAPMl+PAPMl 
APAPM2»APAPM2+PAPM2. 

APPMD1»APPM01*PPMD1 

APPHD2-APPM02+PPM02 

APMLC1»APMLC1*PMLCI 

APMLC2-APMLC2+PMLG2 

AP1»API*P1 

AP2»AP2^P2 

TO CALCULATE S.E. OF AVERAGE ESTIMATORS (WANT S.E. SMALL RELATIVE 
TO DIFFERENCE BETWEEN AVERAGE BIAS! 

XAPMLEl»XAPMLEl»PMLEl*PMLEl 

XAPMLE2-XAPHLE2+PMLE2PPMLE2 

.XAPEPM1»XAPEPMI+PEPMI*PEPM1 

XAPEPM2*XAPEPM2*PEPH2*PEPM2 

XAPAPMI«XAPAPM1>PAPM1*PAPM1 

XAPAPM2«XAPAPH2^PAPM2PPAPM2 

XAPPMOl-XAPPMOl+PPMOlPPPMOl 


ooo r»r»r» ooo ooo nnn oor» r»r»r» ooooo 


XAPPn02«XAPPH02-»PPH02*PPN02 

XAPMLC1«XAPMLCI+PMLC1*PMLC1 

XAPMLC2-XAPMLC2^PMLC2*PMLC2 

XAP1«XAP1+P1*P1 

XAP2«XAP2*P2*P2 

DIFFERENCES OF ESTINATORS FROM EXACT POSTERIOR MEAN 

BIAS OF T.S. APPROX POSTERIOR HEAN FROM EXACT POSTERIOR MEAN 

XAPNNl«XAPNNl>EAPNNl 

XAPHN2*XAPHN2«EAPNN2 

XXAPNN«EAPMN1*EAPNN1«EAPHN2*EAPNN2'»EAPMN3*EAPNN3 

HEAN SQUARE ERROR OF T.S. APN FRON EXACT POSTERIOR MEAN 

XAPNNSQ-XAPNNSQ^XXAPNN 

FOR VARIANCE OF MEAN SQUARE ERROR 

XAPMN4>XAPNN4«XXAPMN«XXAPMN 

BIAS OF POSTERIOR MODE FROM EXACT POSTERIOR MEAN 

XPMOl-XPMDl+EPHOl 

XPMD2»XPM02*EPMD2 

XXPHO»EPM01*EPHD1»EPM02*EPMD2*EPH03*EPMD3 

MEAN SQUARE ERROR OF POSTERIOR MODE FROM EXACT POSTERIOR MEAN 
XPMDSQ>XPMOSQ*XXPMD 

fOR VArIanTe OF “mean SQUARE” ERROR 

XPMO^-XPMOA+XXPNDPXXPMD 

BIAS OF MLE FROM EXACT POSTERIOR MEAN (INCOMPLETE DATA) 

XMLlaXMLl^EMLl 

XML2-XHL24EHL2 

XXNL«ENL1PEML1*ENL2PEML2*EML3*ENL3 

MEAN SQUARE ERROR OF IC>0 MLE FRON EXACT POSTERIOR MEAN 
XMLSQ«XMLSQ4-XXML 
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FOR VARIANCE OF NSE 
XNL4>XNL4«XXNL*XXNL 

DIFFERENCE OF APRRQXIMATED POSTERIOR COVARIANCES FROM EXACT 
POSTERIOR COVARIANCES 

IF (COVSKIP-OI 51»51»50 

50 NCOV-NCOV-1 
COVSKIP-0 
60 TO S3 

CUMULATE VALUES TO AVERAGE EXACT POSTERIOR VARIANCES 

51 AEPMCll-AEPHCll+EPMCll 
AEPMC12-AEPMC12+EPMC12 
AEPHC22>AEPMC22«EPMC22 

CUMULATE VALUES TO AVERAGE APPROXIMATE POSTERIOR VARIANCES 

AAPMCll-AAPMCll+APMCll 

AAPMC12»AAPMC12»APMC12 

AAPMC22*AAPHC22>APMC22 

FOR S.E. OF COVARIANCE AVERAGES 

A2EPC11-A2EPC11*EPMC11*EPHC11 

A2EPC12-A2EPC12+EPMCI2*EPMC12 

A2EPC22-A2EPC22+EPMC22*EPMC22 

A2APC11»A2APC11+^APHC11*APMC11 

A2APC12-A2APC12*APMC12*APMC12 

A2APC22»A2APC22*APMC22*APMC22 

CUMULATE DIFFERENCES BETWEEN EXACT AND APPROXIMATE POSTERIOR 
VARIANCES 

BAPMCll-BAPMCll^EAPMCll 

BAPMC12-BAPMC12+EAPMC12 

BAPMC22«BAPMC22*EAPMC22 

FOR AVERAGE PERCENT RELATIVE DIFFERENCE 

PRDC11»ABS(EAPMC11)/EPMC11 

PRDC12-ABS(EAPMC12/EPMC12) 
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PRDC22«ABS(EAPHC22I/EPHC22 

APRDCll-APROCll^PROCll 

APR0C12*APR0C12«PR0C12 

APR0C22-APRDC22«PR0C22 

FOR nSE OF ELEMENTS OF APPROXIMATED POSTERIOR COVARIANCE MATRIX 

EAP2C11«EAPMC11«EAPMC11 
EAP2C12-EAPMC12*EAPMC12 
EAP2C22>EAPMC22*EAPMC22 
CllMSE* C11MSE^EAP2C11 
C12MSE* C12MSE»EAP2C12 
C22MSE* C22MSE+EAP2C22 
VVV-EAP2C11+EAP2C12+EAP2C22 
VARMSE-VARMSE^VVV 

FOR VARIANCE OF MSE 

VCIINSE-VC11MSE*EAP2C11*EAP2C11 

VC12MSE-VC12MSE*EAP2C12*£AP2C12 

VC22MSE"VC22MSE»EAP2C22*EAP2C22 

VAR4»VAR^4VVV*VVV 

ADO TO BIAS-SIGN AND X-RELATIVE-OIFFERENCE COUNTS 

CALL COUNTStEAPMCll»PROCll»l) 

CALL C0UNTS(EAPMC12»PR0C12>2) 

CALL COUNTS(EAPNC22»PROC22»3» 

ALSO CHECK DIRECTION OF BIAS OF DEPENDENT COVARIANCES 

__EC33-EPMC11 + EPMC22*2_^*^MC12 

EC13— EPMC11-EPMC12 

EC23 — EPMC22-EPMC12 

AC33-APMC11+APMC22*2,*APMC12 

AC13—APMC11-APMC12 

AC23—APMC22-APMC12 

BC33-AC33-EC33 

BC13-AC13-EC13 

BC23-AC23-EC23 

PRDC13«ABS(BC13/ECX3) 

PRDC23«ABS(BC23/EC23) 

PRDC33«ABS(BC33)/EC33 
CALL C0UNTS(BC13*PR0C13t6) 

CALL COUNTS(6C23»PRDC23»7) 
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OF POOR QUALTIYi 


CALL C0UNTS(BC33*PR0C33*8) 


DIFFERENCE OF ESTIMATORS FROM GENERATED P 


93 YMLC1>YMLC1«0MLC1 
YMLC2>YMLC2«0MLC2 

YYMLCO-OMLC1*OMLC140MLC2*OMLC2«OMLC3*DMLC3 

YMLC20-YMLC204YYMLCO 

SQNLCO«YYMLCO*YYMLCO 

REGMLCO-REGMLCO^SQMLCO 

NOTE THAT FOR REMAINING COMPARISONS WITH GENERATED P WE ARE 
MAINLY INTERESTED IN MSE BECAUSE MAIN CONCERN IS DETERMINING WHICH 
ESTIMATOR BEST MINIMIZES QUADRATIC LOSS. HOWEVER* WE WILL ALSO 
CALCULATE BIAS PROM GENERATED lOR GIVEN) P. 

BIAS OF EXACT POSTERIOR MEAN FROM GENERATED P- 

YEPMNl-YEPMNl+DEPMNl 

YEPMN2«YEPMN2*OEPMN2 

YYEPNN«DEPMN1*0EPMN1*DEPNN2*DEPMN2*DEPMN3*DEPNN3 
FOR USUAL MEAN SQUARE ERROR OF EPM FROM GENERATED P 
YEPMSQ-YEPMSQ+YYEPMN 

FOR REGRESSION ESTIMATION AND CONTROL VARIATE MSE EPM 

REGEPMO-REGEPMO+YYEPMNAYYMLCD 

FOR VARIANCE OF ALL 3 MEAN SQUARE ERRORS 

YEPMN4»YEPMNA*YYEPMN*YYEPMN 

BIAS OF APPROX POSTERIOR MEAN FROM GENERATED P 

YAPMNl-YAPMNl+OAPMNl 

YAPMN2-YAPMN2+DAPMN2 

YYAPMN-DAPMNlADAPMNltDAPMN2*DAPMN2«DAPMN3*DAPMN3 
USUAL MEAN SQUARE ERROR FOR APMN FROM GENERATED P 
YAPMSO-YAPMSO+YYAPMN 



ooo ooo nnfi oor» 


C 

C 

c 


FOR RE6RESSI0H ESTIMATION AND CONTROL VARIATE HSE 

REGAPN0*RE6APM0«YYAPMN«YYMLC0 
C 

C FOR VARIANCE OF MEAN SQUARE ERRORS 
C 

YAPMN4»YAPMNA>YYAPHN*YYAPMN 

FOR BIAS OF POSTERIOR MODE FROM GENERATED P 

YPM01«YPH01»0PM01 
YPMD2-YPM02+0PM02 

YYPMO«OPM01«OPM014^0PM02*DPN02-»OPH03*DPMD3 

ypmosq-ypmdso^yypmo 

REGPMOO«REGPMDO^YYPMO*YYMLCO 
YPMOA»YPMO^^YYPMO*YYPMO 

BIAS OF MLE FROM GENERATED P 

YHll-YMLl+OMLl 
YML2-YML2+OML2 

YYML«OMU*OHL1*OML2*OHL2^0HL3*OML3 
YNLSQ»YMLSQ>YYML 
REGNLO«REGNLO+YYHl*YYMLCO 
YHL4-YMLA+YYML*YYML 

BY TWO OIFF CRITERION CALC BEST EST FOR THIS TRINOMIAL TRIAL 
CALL BESTEST(BESTQL(1)> 



WAPMN1«WAPMN1«EAPMN1*EAPMN1 

WAPMN2»WAPHN2*EAPHN2*EAPMN2 

WPMD1«WPM01*EPMD1*EPM01 

MPM02»WPM02*EPMD2*EPMD2 

WML1>WML1«EML1«EML1 

WML2«WML2>EML2«EML2 

OMLCl»UMLCH-OMLCl*OMLCl 

UMLC2«UMLC2*0MLC2*0HLC2 

UEPHNl-UEPMNl+OEPMNlPOEPMNl 

UEPMN2-UEPMN2*0EPMN2*0EPMN2 

UAPMN1«UAPMN140APMN1*0APMN1 

UAPMN2«UAPMN2*DAPMN2*DAPMN2 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


UPM0l«UPM01*0PM01*DPM01 

UPM02-UPN02*0PMD2*DPM02 

UHL1»UML1>0ML1*DML1 

UHL2*UNL2^OHL2«0ML2 

CALCULATE ROBUSTNESS ESTIMATORS (FOR OUADRATIC>LOSS COMPARISON 
ONLYI. USE WRONG PRIOR IN ESTIMATES. 

EXCLUDE EPM BECAUSE OF EXPENSE. 

ROBUSTNESS SET 1. UNIFORM PRIOR. RECALCULATE ONLY APPROXIMATE 
POSTERIOR MEAN BECAUSE FOR A UNIFORM PRIOR THE POSTERIOR MODE 
WILL EQUAL THE ALREADY CALCULATED M.L.E. (FOR INCOMPLETE DATAI 

Z1N>Z14-1. 

Z2N*Z2«1. 

Z3N«Z2>1. 

SSN-SS43. 

IROBUST-1 

CALL METHODS 

IF (ISTOP-1) 55»58»58 

55 VAPMNI*VAPMN1>DAPMN1 
VAPHN2«VAPMN2*DAPMN2 

VVAPMN«0APHN1*DAPMN1>DAPMN2*0APMN2*DAPMN3PDAPMM3 

AVPMN1«AVPMN1*PAPM1 

AVPMN2*AVPHN2+PAPM2 

FAPHN1*FAPMN1+0APMNI*DAPMN1 

FAPMN2«FAPMN2*OAPHN2ADAPMN2 

AFPMNl*AFPMNltPAPMl*PAPMl 

AFPMN2«AFPMN2+PAPM2*PAPM2 ' 

VAPMSQ-VAPMSQ+VVAPMN 

REGAPMI«REGAPM1>VVAPMNPYYMLC0 

VAPMNA-VAPMNA+VVAPMNAVVAPMN 

YMLC21-YHLC21+YYMLCD 

REGMLCI-REGMLCI+SOMLCD 

BY TWO DIFF CRITERION CALC BEST EST FOR THIS TRINOMIAL TRIAL 
CALL BESTEST(BESTQL(9)I 

ROBUSTNESS SET 2. IN BAYESIAN ESTIMATORS APMN AND PMD« USE 
PRIOR PARAMETERS 10. (NU/IO. ♦(0.09*0. 0S*-0.1An 

56 ZlN>ZltXNU1^0.9 
Z2M»Z2+XNU2^0.5 
Z3M»Z3+XNU3-1.4 
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ssn-ss+xnu 

IR0BUST*2 

CALL METHODS 

IF CISTOP-1) 57»58»58 

57 QAPMN1>QAPMN1>0APMN1 
QAPHN2-0APHN240APMN2 
AQPMN1-AOPMN14PAPM1 
A0PHN2"AQPMN24PAPH2 

QOAPNN-DAPMNl«0APMNl40APMN2*DAPnN2*DAPHN3*0APMN3 

QAPMSQ-0APMS04Q0APMN 

REGAPM2»REGAPM24QQAPMN*YYMLC0 

0APnN^«QAPMN4400APMN*Q0APMN 

YMLC22-YMLC224YYMLCD 

REGHLC2-RE6NLC24SQHLC0 

OPN01-QPM0140PM01 

QPM02"QPM02*0PM02 

AQPMDI-A0PM014PPMD1 

AOPMD2-AQPM024PPM02 

QQPMO«OPM01*OPM0140PM02*OPM0240PM03*DPM03 

OPHOSO-QPNOSQ4QOPMO 

REGPM02«REGPM02400PMO*YYMLCO 

OPH04«QPH044QQPHO*QQPHO 

GAPMN1-GAPHN140APHN1*0APMN1 

6APMN2-GAPMN240APMN2P0APMN2 

GPM01-GPM0140PM01*DPMD1 

GPMD2«GPM02*0PM02*0PM02 

AGPMN1«AGPHN14PAPM1*PAPM1 

AGPHN2«AGPHN24PAPM2*PAPM2 

PGH0W6M0l4PPM0l*P-PM01 

PGMD2»PGM024PPM02*PPM02 

BY TWO OIFF CRITERION CALC BEST EST FOR THIS TRINOMIAL TRIAL 
CALL BESTEST(BESTQL(17) I 

GO TO 65 

58 ISTOP-0 

IF (IROBUST'll 60»61>62 
60 NUMXZ-NUMXZ-1 
NC0V»NC0V-1 
NXZRl-NXZRl-1 
NXZR2-NXZR2-1 
GO TO 65 
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original Pacj}' to 

^OOR QU^ 


61 NXZRl-NXZRl-1 
60 TO S6 

62 NXZR2-NXZR2-1 
65 CONTINUE 

AVERAGE OVER NUNXZ TRINOMIAL RESULTS FOR FIXED XNU>(XNU1»XNU2*XNU3 
GENERATED P-(R1»P2»P3)# PI0» AND SS 


ESTIMATOR MEANS (ESTIMATORS AVERAGED OVER NUMXZ TRIALSI 

APMLEl'APMLEl/NUMXZ 

APNLE2-APMLE2/NUMXZ 

APEPMl-APEPMl/NUMXZ 

APEPM2-APEPM2/NUMXZ 

APAPMl-APAPMl/NUMXZ 

APAPM2-APAPM2TNUMXZ 

APPMOl-APPHDl/NUMXZ 

APPN02«APPM02/NUNXZ 

APMLCl-APMLCl/NUHXZ 

APMLC2-APMLC2/NUMXZ 

APl-APl/NUHXZ 

AP2-AP2/NUMXZ 

AVPHNl-AVPMNl/NXZRl 

AVPMN2-AVPMN2/NXZR1 

AQPMN1-A0PHN1/NXZR2 

AQPMN2-AOPMN2ZNXZR2 

AQPM01-AOPM01/NXZR2 

AQPMD2-A0PM02/NXZRZ 

N0«NUMXZ*(NUMXZ-1» 

NC«NC0V*(NC0V-1) 

NDRl-NXZRl*(NXZRl-l) 

NDR2»NXZR2*(NXZR2-1) 

standard errors of ESTIMATOR MEANS 

SEl-SQRTI (XAPNLE1-NUMXZ*APMLE1*APMLEI)/ND) 
SE2*S0RTnXAPNLE2-NUMXZ*APMLE2*APHLE2»/ND» 
SE3«S0RT{(XAPEPM1-NUMXZ*APEPMI*APEP«1)/ND) 

SEA*SORT( (XAPEPMZ-NUMXZ*APEPM2*APEPM2)/ND) 

SE5»S0RT( ( XAPAPM1-NUMXZ*APAPM1PAPAPM1»/N0) 

SE6«S0RT( (XAPAPM2-NUMXZ*APAPM2*APAPM2)/ND) 

SE7«SQRT( (XAPPM01-NUMXZ*APPH0l*APPH01) /ND) 
SE8»SQRT(<XAPPN02-NUMXZ*APPN02*APPMD2»/N0) 
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SE9«SQRTnXAPMLCl-NUMXZ*APHLCl*APHLCll/N0» 

SEIO-SORTC <XAPMLC2-NUMXZ*APMLC2*APMLC2)/N0» 

SE11»SQRT( (XAP1-NUMXZ*AP1*AP1»/N0» 

SE12*SORT( (XAP2-NUNXZ«AP2*AP2)/N0) 

SE13-S0RT< <AFPMN1-NXZR1*AVPMN1*AVPMN11 /N0R1> 

SE14«S0RT( I AFPMN2>NXZRl«AVPHN2*AVPnN2)/NDRl) 

SE15«SORT( ( AGPMN1-NXZR2«AQPHN1«AQPNN1) /NDR2) 

SE16«S0RT( (AGPNN2-NXZR2*AQPHN2*AQPNN2) /NDR2) 

SE17»SQRT( <PGN01-NXZR2*AOPM01*AQPNOU/NOR2) 

SE18«SQRT( (PGMD2>NXZR2PAQPH02PAQPH02)/N0R2) 

PRINT 1079* NUMXZ 

PRINT 1080* APNLE1«SE1>APEPN1>SE3*APAPN1»SE5»APPHD1>SE7»APNLC1»SE9 
l*APl«SEllf APHLE2*SE2«APEPN2«SE^*APAPH2»SE6*APPN02»SEa#APNLC2»SE10# 
2AP2*SE12 

1079 FORHAT<//* AVERAGE ESTIMATORS (MEANS) AND THEIR STANDARD ERRORS OV 
1ER**IA»* TRIALS*//) 

1080 F0RMAT(10Xf*APMLE**7X»*S.E«**6X**APEPM*»7X»*S.E.*»6X»*APAPM*»7X»PS 
l.E.*»6X»*APPM0*»7X»*S.E.*»6X»*APMLC*»7X»*S.E.*#9X»*AP*»7X»*S.E.*/* 
2 PI *12F11.7/* P2 ♦12F11.7/) 

AVERAGE BIASES AND COVARIANCE ESTIMATORS 

XAPMNl-XAPMNl/NUMXZ 
XAPMN2-XAPMN2/NUMXZ 
XPM01«XPM0I/NUMXZ 
XPM02»XPM02/NUMXZ 
XMLl-XMLl/NUMXZ 
XML2-XML2/NUMXZ 
. AEPMC11«AEPMC11/NC0V 
AEPMC12-AEPMC12/NC0V 

AEPNC22-AEPMC22/NC0V 

AAPMCll-AAPMCir/NCOV 

AAPMC12-AAPHC12/NC0V 
AAPMC22-AAPMC22/NC0V 
. BAPMCll-BAPMCll/NCOV 
BAPMC12-BAPMC12/NCOV 
BAPHC22-8APMC22/NC0V 
YMLCl-YMLCl/NUMXZ 
YMLC2-YMLC2/NUMXZ 
YEPMNl-YEPMNl/NUMXZ 
YEPMN2-YEPMN2/NUMXZ 
YAPMNl-YAPMNl/NUMXZ 
YAPMN2-YAPMN2/NUMXZ 
. YPMOl-YPMOl/NUMXZ 
. YPM02-YPM02/NUMXZ 
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YMLl-YMLl/NUMXZ 

YML2-YML2/NUMXZ 

VAPHMl-VAPHNl/NXZRl 

VAPHN2-VAPHN2/NXZR1 

OAPNN1-OAPHN1/NXZR2 

OAPHN2-QAPHN2/NXZR2 

QPN01«0PMD1/NXZR2 

QPN02-QPH02/NXZR2 

EBIASKNREPLIC* ISS»IPI0I«XAPHN1 

EBIAS2(NREPLIC*ISS*IPI0)-XPH01 

EBIAS3(NREPLIC*ISS«IPI0)>XHL1 


CALCULATE S.E. OF BIAS. WANT THIS TO BE SMALL RELATIVE TO THE 
DIFFERENCE BETWEEN THE BIASES 


SE21-SORT( 

SE22«SQRT( 

SE23»SORT( 

SE24«SQRT( 

SE25-SQRTI 

SE26-SQRT( 

SE27-SORT( 

.SE28-SQRT( 

SE29-SQRT( 

SE30-S0RT( 

SE3l«SORT( 

SE32-SQRT( 

SE33«S0RT( 

SE3A«SORT( 

SE35«SQRT( 

SE36>SORT< 

SE37»SQRT( 

SE38«SQRT( 

SE39«S0RT( 

SEAO>SQRT< 

SE41-SQRT( 

SE42»SaRT( 

$E43«SQRT( 

SE44»SQRT( 

SE45»SQRT( 

SE46-SQRT( 

SE47»SQRT( 

SE48>SORT( 

SE49«SQRT( 

SE90-SQRT( 


(WAPMN1>NUMXZ*XAPHN1*XAPMN1) /NO) 

(WAPMN2-NUMXZ*XAPMN2«XAPMN2)/ND> 

<WPH01-NUMXZ«XPMDl*XPM0l) /ND> 

(WPM02-NUMXZ*XPM02*XPM02> /NO) 

(WMLl-NUMXZ*XMLl*XMLl)/ND» 

(WML2-NUMXZ*XML2*XML2)/N0) 

{A2EPCll-NCOV*AEPMCll*AEPMCll)/NC» 

<A2EPC12-NCOV*AEPMC12*AEPMCI2)/NC» 

(A2EPC22-NCOV*AEPMC22*AEPMC22» /NO 

(A2APCll-NCOV*AAPMCll*AAPMCll ) /NO 

(A2APC12-NC0V*AAPHC12*AAPHC12 )/NO 

(A2APC22-NCOV*AAPHC22*AAPMC22) /NO 

(C11MSE-NC0VPBAPMC11*8APMC11> /‘NO 

<C12HSE-NC0V*BAPMC12*BAPMC12 )/NC ) 

(C22MSE-NC0V*BAPHC2Z*8APMC22) /NC » 

<Uf1LCl-NUMXZ*YMLCl*YHLCl) /NO) 

(UMLC2-NUHXZ*YMLC2*YMLC2»/N0) 

(UEPMN1-NUMXZ*YEPMN1*YEPMN1) /NO) 

(UEPMN2-NUNXZ*YEPMN2*YEPMN2)/ND) 

(UAPMNl-NUMXZ*YAPHNl*YAPMNl)/ND) 

(UAPMN2-NUMXZ*YAPMN2*YAPMN2)/ND) 

(UPM01-NUHXZ*YPM0l*YPMDl)/ND) 

(UPM02-NUMXZ*YPMD2*YPM02»/N0) 

<UHL1-NUHXZ*YML1*YM11) /NO) 

(UML2-NUHXZ*YHL2*YML2)/ND> 

(FAPMN1-NXZRI*VAPHN1*VAPMN1)/NDR1> 

CFAPMN2-NXZR1*VAPHN2*VAPHN2)/NDR1» 

<GAPMNl-NXZR2*QAPMNl*0APNNl) /N0R2» 

(GAPMN2-NXZR2*0APMN2*0APMN2)/NDR2> 

(GPMD1>NXZR2*QPMD1*QPM01)/NDR2) 
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SES1*SQRT< (GPMD2-NXZR2*QPM02*QPN02}/NDR2) 

C 

PRINT 1150* NUMXZ*XAPMNl*SE21*XPNDl*SE23*XNLl*SE25*XAPHN2*SE22*XPn 
lD2»SE2^»XML2f SE26*yEPMNl»SE30»YAPNNl,SEAO*yPMDl,SEA2»YMLl*SEA<*»YML 
2Cl»SE36f YEPMN2»SE39*YAPMN2»SE41»YPM02»SEA3»YML2»SE45»YMLC2*SE37 

1150 FORHATI//* AVERAGE BIASES AND THEIR STANDARD ERRORS OVER**IA»* TRI 
lALS*//* FROM EPM*15X*XAPM*16X*S.E.*16X*XPMD*16X*S.E.*17X*XML*16X*S 
2*E.4-/10X**P1 ♦»6F20.1A/10X»*P2 «* 6F20.1A / / * FROM P*<»X» *YEPM*8X*S . 
3E.*8X*YAPM*8X*S .E.*8X*YPMD*0X*S.E.*9X*YML*8X*S.E.^8X*YMLC*8X*S.E.* 
A/10X*P1 ♦10F12.10/10X+P2 *10F12.10/) 

PRINT 1154* AVPMN1*SE13» APMLE1*SE1»AQPMN1>SE15*A0PM01>SE17*AVPNN2* 
1SE1A»APMLE2»SE2.AQPMN2,SE16.AQPM02»SE18 

1154 FORMAT!//* RQBUST'ESTIMATOR AVERAGES AND S.E.^S. WANT S.E. SMALL 
IRELATIVE TO DIFFERENCE BETWEEN ESTIMATORS .♦/ 9X •APMR1*11X *S . E . *9X*P 
2MDR1»MLE*8X*S.E.*11X*APMR2*10X*S.E.*11X*PMDR2*10X*S.E.*/* PI *4! FI 
32.8»F18.14»/* P2 *4 ( F12. 8»F18 .14 ) // ) 

PRINT 1155* VAPMN1*SE46* YML1»SE44*0APMN1*SE48»QPND1*SE50*VAPMN2»SE 
147* YML2*SE45*0APMN2*S£49,0PM02*SE51 

1155 FORMAT!//* ROBUST-ESTIMATORS BIASES AND STANDARD ERROR S . */10X* * 

1VAPM*11X**S.E.**9X**VPM0-YML*9X**S.E.*11X*6aPM*11X*S.E.*11X*0PM0*1 
21X*S.E.*/* PI *4!F12.8*F18.14J/* P2 *4 ! F12 .8* F 18 . 14) // ) 

MEAN SQUARE ERRORS AND THEIR STANDARD ERRORS 

XAPMNSQ-XAPMNSO/NUMXZ 
XPMDSO-XPMDSO/NUMXZ 
XMLSQ-XMLSO/NUMXZ 
EMS1!NREPLIC.ISS*IPID)-XAPMNSQ 
EHS2!NREPLIC*ISS* IPIO)-XPMOSQ 
EMS3!NREPLIC*ISS*IPI0)»XMLS0 
YMLC20-YMLC20/NUMXZ 
YEPMSQ-YEPMSO/NUMXZ 
' Y”APHS'0«YA>M‘50/NUMXZ 
YPMOSQ-YPMDSO/NUMXZ 
YMLSQ-YMLSQ/NUMXZ 

SDl-SORT! !XAPMN4-NUMXZ*XAPMNSQ*XAPMNS0)/ND) 

SD2-S0RT! ! XPMD4-NUMXZ*XPMDSQ*XPMDSQ) /ND) 

SD3-S0RT! ! XML4-NUMX Z*XML SO*XHLSO ) /NO ) 

504- SQRT! !REGMLC0-NUMXZ*YMLC20*YMLC20) /NO) 

505- SQRT!! YEPMN4-NUMXZ*YEPMSQ*YEPMS0) /NO) 

506- SQRT! !yAPMN4-NUMXZ*YAPMS0*YAPMS0)/ND) 

507- S0RT!! YPMD4-NUMXZ*YPMOSO*YPMDSQ) /NO) 

508- SQRT!! YML4-NUMXZ*YMLSQ*YMLSQ) /NO) 

C 

C DIFFERENCES BETWEEN MEAN SQUARE ERRORS 


original pagp tq 

0^ POOB 


OHl-XAPMNSO-XPMDSQ 

0N2-XAPMNSQ-XHLS0 

0«3»XPH0Sa-XMlS0 

0M^«YMLC20-YEPHSQ 

0M5"YfaC20-YAPMSQ 

0H6-YHLC20-YPM0S0 

0H7-YHLC20-YMLSQ 

DMB-YEPMSO-YAPMSO 

0M9-YEPMS0-YPM0SQ 

OMIO-YEPMSO-YMLSQ 

OHll-YAPMSQ-YPMOSO 

DM12-YAPMSQ-YHLS0 

OM13-YPMOSO-YMLSO 

PRINT 1200a NUHXZ»XAPMNSO»S01»YEPHSQ*SD5aYMLC20*SD4»XPMDS0*SD2»YAP 
IHSQ»S06aYMLS 0»S08»XMLSO*S03»YPMDSO»S07 
1200 FORMAT!//* AVERAGE MEAN SQUARE ERRORS AND THEIR STANDARD ERRORS OV 
1ER*IA* TRIAL S*//17X,*MSE*12X*S.E.*25X*MSE*12X*S.E.*25X*MSE*12X*S.E 
2.*/2X*APM-EPM*F14.9»F18.13»7X»*EPM-P*FlA,9»F18.13»7X»*MLC-P*FlA.9» 
3F18.13/2X*PM0-EPM*FIA.9»F18.13»7 Xa*APM-P*F1A.9aF 18.13»8X»*ML-P*FIA 
A.9»Fl8.13/3Xf ♦ML-EPM*F1A.9»F18.13 a7X,*PMD-P*F1<..9»F18.13/) 

PRINT 1300» OHI»OM8»DM11»OM5»OM2aDM9,DM12aOM 6»DM3»DMIO»DM13»DM7 
1300 FORMAT!//* DIFFERENCE BETWEEN AVERAGE MSE"S. WANT DIFFERENCE LARG 
XE RELATIVE TO SEIMSEJ*//* EAPM-EPM0*F1A.9»12X*DEPM-0APM*FIA.9»12X 
2*0APM-0PM0*F1A.9»12X*DMLC-0APM*F19,9/* EAPM-ENLE*F1A.9» 12X*DEPM-D 
3PM0*F1A«9»12X*0APM-0MLE*F19.9»12X*0MLC-DPMD*FIA.9/* EPM0-EMLE*F1A 
9.9»12X*OEPH-OMLE*F19,9»12X*OPMD-DMLE*FIA,9»12X*DMLC-DHLE*F1A.9/) 
XN-NUMXZ 
DO 71 I-l»6 
IF !I-A) 69»67»68 

67 XN-NXZRl 
60 TO 69 

68 XN-NXZR2 

69 AVNUMIT!Il»AVNUMIT!l l/XN 
DO 70 J>1 a10 

70 CTNUMIT!I» J)-CTNUMIT!If J )/XN 

71 CONTINUE 

PRINT ZOOOf IPaPI0aSS»NREPLICaNUMXZ*NXZR1*NXZR2*!AVNUMIT!I)»I>I*6) 
2000 FORMAT!///* NUMBER OF ITER FOR CONVERGENCE AVERAGED OVER NUMBER OF 
I TRINOMIAL SIMULATIONS. IP»*I3»* PI0-*FA,2»* SS-*F3.0»* NREPLIC-* 
212/* NUMXZ-*I3»* NXZR1«*I3»* NXZR2"*I3»* AV NUN ITER FOR MLE«*F7 
3.3* PM0R0«*F7.3* APMR0«*F7.3* APMR1»*F7.3* PMDR2-*F7.3* APMR2«*F7. 
A3//I 

PRINT 2010a !!CTNUMIT!Ia J lAJ*lAlO)Al-l»AA3)A!!CTNUMIT!I»J)«J«lrlO) 
1aI-2a5#3)a !!CTNUM1TIIa J) aJ-1a10)a1-3a6a3) 



o o o o 


2010 F0RHAt(« PROPORTION OF DATA SETS FOR WHICH NUMBER OF ITERATIONS WA 
IS OF SPECIFIED AM0UNTS*/10XP1 23 A 5 6 7 

2 8-10 11-15 GT 15*11X*1 2 3 A 5 6 7 B 

3-10 11-15 GT 15*/* MLE *10F6.3* APMRl *10F6.3/* PMDRO *10F6,3* 

A PHDR2 A10F6.3/* APMRO A10F6.3* APMR2 *10F6.3//J 
YEPMSQ»NUMXr*YEPMSQ 
YAPMSO«NUHXZ*YAPMSQ 
YPHDSO-NUMXZ*YPMDSO 
YMLSQ-NUMXZPYHLSO 
. YMLC20-NUMXZ*YMLC20 

T*C1.-IP1*P1+P2*P2*P3*P3)>/SS 

C 

CALL ESTMSE<YEPMSQ»Y£PMNA»REGEPH0»YMLC20»REGMLC0»T»NUMXZ»MSE(1) > 
CALL ESTMSE(YAPMSO» YAPMNA»REGAPM0»YMLC20»REGHLC0» T»NUMXZ»MSE (7> » 
CALL ESTMSE(YPM0SQ»YPM0A»REGPM00»YMLC20»REGMLC0»T»NUMXZ»MSE{13>I 
CALL ESTMSE(YMLSQ»YMLA*REGML0*YMLC20»REGNLC0»T»NUMXZ,HSE(1<»») 

CALL ESTMSEtVAPMSQ# VAPMNA»REGAPMl,YMLC21»REGMLCl»T»NXZRI»MSE<25n 
CALL ESTMSE(0APMS0»0APMNA,REGAPM2»YMLC21»REGMLC2»T>NXZR2»HSE<31> > 
CALL ESTMSE(QPM0SQ»QPM0A»REGPM02»YMLC22#RE6MIC2»T»NXZR2»HSE(37)) 

C 

PRINT 2030* IPtPIOf SS»NREPLIC» ((MSE(I»J)»J*1»7)»I>1»6) 

2030 FORMAT!* THREE KINDS OF MSE (AND THEIR VARIANCES! FOR OUADRATIC-LO 
ISS COMPARISONS. IP«*I2f* P1D«*FA.2,* SS«*F3.0»* NREPLIC-*I2//18X* 
2EPM*15X*APM*15X*PM0*15X*MLE*IAX*APMR1*13X*APMR2*13X*PMDR2*/* REG M 
3SE *7E18,7/* VAR(MSEI*7E10.7/* CV MSE ♦7E18.77* VAR < MSE » *7E 18. 77* 

A RE MSE ♦7E18.77* VAR ( MSE) *7E18.777 > ^ 

C 

IF (ISS-1) 2035>2035.20AO 
2035 0LMS11(IGEN>1PID.NREPLIC)«MSE(5»2) 
OLMS21(IGEN»IPIO»NREPLIC)-M$E(5>3) 
0LMS3UI6EN»IPI0»NREPLICI-MSE(5.A) 

OL M S A1 ( I GE N # I P I D » N REP L I C ) • M S E (5,5 > 

0LMS51IIGEN,IPI0,NREPLIC)»MSE(5,6) 

0LMS61(IGEN,IPI0,NREPLIC)-MSE(5,7) 

GO TO 20A5 

20A0 0LMS12(IGEN,IPID,NREPLIC)«M$E(5,2) 
QLMS22(1GEN,1PID,NREPL1C)«MSE(5,3) 

0LMS32(IGEN, IPID,NREPLIC)-MSE(5,A) 

0LMSA2(IGEN,IPI0,NREPLIC) -MSEC 5,5) 
QLMS52n6EN,lPI0,NREPLIC)«MSE(5,6) 
OLMS62(I6EN,IPIO,NREPL1C)«MSE(5,7) 

PROPORTIONS FOR BEST ESTIMATOR (BEST IN TERMS OF SMALLEST SUMMED 
SQUARED ERROR AND X REL DIFF FOR SUM BEING OVER THE THREE 
COMPONENTS OF AN ESTIMATOR) AND FOR SIGN OF BIAS 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


c 

2045 CTSEOLCD-NUNXZ 
CTSEQL(2)«NXZK1 
CTSEQL(3)-NXZR2 
DO 75 I-l»4 
DO 72 If(«l*3 

BESTQL(I»1« IR)-BESTQL(I*l>IR)/CTSEOl(IR) 

BESTQL(I«2» IR)«eESTQL(I*2*IR)/CTRD0LIIR) 
SBIASQL(IR»I)-S8IAS0L(IR*n/NUNXZ 

72 CONTINUE 
IF II-4) 73#75»75 

73 BESTEF(I*1)-BESTEP(I>1)/NUHXZ 
BESTEPCI»2»-BESTEP{I*2I/CTR0EP 
DO 74 K»l»3 

SBIASEP(K>n*SBIASEP(KfI I/NUNXZ 

74 CONTINUE 

75 CONTINUE 
DO 76 K»l»3 

SBIASQL(K»5)*SBIAS0L(K«5)/NXZR1 
SBIAS0L(K,6»»SBIAS0L(K»6WNXZR2 
S8IASQL(K*7I-SBIASQL(K>7I/NXZR2 

76 CONTINUE 

CALCULATE Z ABS REL OIFF LESS THAN (INSTEAD OF BETWEEN) SPECIFIED 
AMOUNTS 

DO 80 I-l»7 
DO 79 II«2»0 

PROQL(II»I)-PROQL(II>I)>PROQL(II>1>I) 

IF (1-4) 78 , 79,79 

78 PROEP(II*I)-PROEP(II»I)*PROEP(II-l*I) 

79 CONTINUE 
BO CONTINUE 

IR-l 

DO 85 I>1»7 
IF (1-5) 83«81#82 

81 IR>2 
60 TO 83 

82 IR>3 

83 DO 84 II-I>8 

PRDQL(II*I)-PRDQL(II»I)/CTRDQL(IR) 

84 CONTINUE 

85 CONTINUE 
DO 88 I«l»3 
DO 88 II>1>8 



r» r» r» r> o o r» r> o 


26 


PRDEPdl*! l•PROEP(II»n/CTROEP 
86 CONTINUE 

PRINT 2050* ((BESTEP(I*J)tI«I»3)* flBESTOLd* J*K|*I«2*4)*K>1*3)* J> 
11 * 2 ) 

2050 FORHATt* PROPORTION OF CASES THAT AN ESTIMATOR IS BEST. FIRST 3 C 
lOLNS ARE RESULTS FOR ESTIMATING EPM. REMAINING COINS* FOR MIN QUA 

20 LOSS. «/A3X«0RIG. -PRIOR ROBUST SET UNIFORM-PRIOR ROBUST SET 

3 PERTURfl-PRIOR ROBUST SET*/20X*MLE PMD APM*11X* *MLE PM 
AD APM*1AX*MLE • PMO APM*13X*MLE PMD APM*/* SOD ERR CRIT 
5 «3F6.2*8X3F6.2*3X*2(9X«3F6.2)/* REL DIFF CRIT «3F6. 2* BX» 3F6. 2* 3 
6X*2(9X*3F6.2)//) 

PRINT 2070* ((SBIASEPd«J)*J«l*3)*(SBIASQL(I*K)*K»l*7)*I-l*3) 

2070 FORMAT!* PROPORTION OF CASES IN WHICH DIFFERENCE BETWEEN FIRST COM 
IPONENT OF ESTIMATOR AND THAT OF ESTIMATED IS OF A CERTAIN SIGN*//1 
25X*EMLE*6X*EPMD*6X*EAPM*16X*QEPM*6X*QMLE*5X*0PMDR0*AX*0APMR0*9X*0A 
3PMRl*9X*0PM0R2*AX*QAPMR2*/6X*NEG *3F10. A* lOX* AFIO. A* 5 X* F 10. A* 5X* 2F 
AIO.A/5XAZERO *3F10.A*10X*AF10,A*5X»F10.A»5X*2F10.A/6X*P0S ♦3F10.A* 
S10X*AFIO.A*5X*F10.A*5X»2F10.A//) 

PRINT 2060* ((PR0EPdI*II*I«l»3)*(PR00LdI*K)*R>l»7)*II«l*8) 

2080 FORMAT!* PROPORTION OF CASES FOR WHICH Z ABSOLUTE RELATIVE DIFFERE 
INCE FOR EACH OF THE THREE ESTIMATOR COMPONENTS IS LESS THAN SPECIF 

21 ED AMaUNTS*/15X*EMLE*6X*EPM0*6X*EAPM*I6X*QEPM*6X*QMLE*5X*QPMDR0*A 
3X*QAPMR0*9X*0APMR1*9X*0PMDR2*AX*QAPMR2*/6X*0.01*3F10.A»10X»AF10.A* 
A5X*F10.A*5X*2F10.A/6X*0.1 *3F10. A* lOX* AF 10. A* 5X* F 10 . A* 5X * 2F10 , A/6X 
5*1.0 *3F10. A*10X*AF10.A»5X* F10.A*5X»2F10,A/5X** 5.0 *3F10.A» 10X*AF 
610.A*5X*F10.A*5X»2F10.A/5X**10.0 *3F10. A* lOX * AF 10 . A* 5X * FIO. A* 5X> 2F 
710. A/ 5X* *15.0 ♦3F10.A,10X,AF10.A*5X*F10.A*5X»2F10.A!5X*20.0 *3F10. 
8A*10X*AF10.A*5X*F10.A*5X*2F10.A/5X**25.0 *3F10. A* lOX* AFIO. A* 5X*F10 
9«A*5X*2F10.A/7) 

PERCENT AVERAGE RELATIVE DIFFERENCE FOR COVARIANCE ESTIMATES 

PAR0C11«100.*!AEPMC11-AAPMC11I/AEPMC11 
PAROC12-100.*(AEPMC12-AAPMC12)/AEPHC12 
PAR0C22-10O.*!AEPMC22-AAPMC22)/AEPMC22 

AVERAGE PERCENT RELATIVE DIFFERENCE 

APR0C11>100.*APR0C11/'NC0V 
APRDC12»100.*APR0C12/NC0V 
APRDC22«100.*APR0C22/NC0V 

SQUARE ROOT MSE DIVIDED BY AVERAGE EPV 

CllMSE-CllMSE/NCOV 
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C12MSE-C12MSE/NC0V 
C22MSE-C22MSE/NCQV 
CllRTl-S0RT(Cllf1SE)/AEPMCll 
C12RT1«SQRT(C12MSE>/AEP«C12 
C22RT1«S0RT(C22«S£)/AEPMC22 

SE(HSEI/HSE 

SE11HSE«S0RT<(VC11MSE- NC0V*C11«SE*C11MSEI/NCI 
SE12MSE»S0RT((VC12MSE- NC0V*C12MSE*C 12MS E ) /NC ) 

SE22MSE«SQRT< (VC22MSE- NC0VPC22rtSE*C22NSE ) /NC ) 

CllRT2-SEllMSE/ClinSE 
C12RT2-SE12MSE/C12MSE 
C22RT2-SE22NSE/C22HSE 
00 91 J*l»8 
00 90 I«l*6 

COUNTROUf J)«C0UNTRDn»J)/NC0V 
IF <I-3» 89»fl9»90 
69 COUNTen*J)*COUNTB<I*J)/NCOV 

90 CONTINUE 

91 CONTINUE 
C 

PRINT 3000# AEPMC11»SE27»AAPMC11#SE30»PARDC11»APR0C11»C11RT1#C11RT 
12»AEPNC12*SE28f AAPNC12»SE31»PARDC12»APRDC12 fC 12RT1>CI2RT2>AEPNC22» 
2SE29»AAPNC22»SE32»PAROC22»APR0C22»C22RTl»C22RT2 
3000 F0RMAT(////12X**AVERAGE EPV*8X*S.E.*9X*AVERAGE APV*8X*S,E.*7X*T. AV 

1 REL OIFF AV Z REL OIFF SQRTIHSEJ /(I » SE (NSE » /NSE*//3X*Cll . 
2P8E16,7/3X*C12 *8E16.7/3X*C22 *8E16.7» 

PRINT 3005» NC0V»IP»PI0,SS»NREPLIC»(C0UNTRD(I»l)»I«l»6)» ICOUNTROt 
1I*3)*I*1*6)»(C0UNTR0(I»8)>I-1>6)»(CDUNTRD(I>2)»I-1»6)»(C0UNTR0(I»6 
2)*I-l*6l*(CQUNTR0(If 7)»I«1»6) 

3005 FORMAT(///P PROPORTION OF*IA* CASES IN WHICH PERCENT REL OIFF WAS 
UESS THAN VARYING AMOUNTS. IP»*I2* PID-*FA.2* SS"*FA.O* NREPLI 
2C-*I2//8X*.01 O.l 1,0 5.0 10. 15,*12XP.01 0.1 1.0 

35.0 10. 15.*12X*.01 0.1 1.0 5.0 10. l5.«/« C11P6F6. 

A3»6X*C22P6F6.3»6X*C33*6F6.3/* C 12*6F6. 3» 6X*C 1 3*6F6, 3/ 6X*C23*6F6. 3 
5//) 

PRINT 3010» NC0V»(C0UNTB(I»ll*I*l>3}»fC0UNTB(I»3)»I>l»3)*(C0UNTB(I 
1>8I»I>1>3)» (COUNTS ( I»2 ) f I-1>3I» (C0UNTBfI»6)» I-1#3)*(C0UNTB( I»7I»1- 
21>3) 

3010 FORMAT!* PROPORTN OF*IA* CASES IN WHICH BIAS IS OF A CERTAIN SIGN* 
1//7X*NEG ZERO P0S*7X*NEG ZERO P0S*7X*NEG ZERO P0S*7X*NEG 

2 ZERO P0S*7X*NEG ZERO P0S*7X*NEG ZERO POS*/* C11*3F6.3* C2 
32*3F6.3* C33*3F6.3* C12*3F6.3* C13*3F6,3* C23*3F6.3I 
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C AVERAGE TRUE PERCENT (PROPORTiaN) INCQNPLETE DATA 
C 

AVTPID-AVTP104TPIO 

AVOPID*AVOPIO«OPID 

PRINT 3070* PIO*NUNXZ»AVTPIO#AVOPIO 

3070 FORNATIP GIVEN PIO IS*FA.2P AVERAGE OVER NUNXZ«*I3P TRIALS OF TRUE 

1 PIO ISPF6.2* AV OIFF BETWEEN TRUE AND GIVEN PID OVER THESE TRIALS 

2 ISPF6.2/////) 

C 

95 CONTINUE 
C 

KASE>KASEPl 

60 TO (9902*9903»9904*9905) KASE 

9902 PIO-0.40 
IPlb«2 
IPRINT-l 
60 TO 1 

9903 PIO-0.15 
NSS-50 
IPIO-1 
ISS-2 

GO TO 1 

9904 PIO-0.40 
IPID-2 
ISS-2 

GO TO 1 

9905 CONTINUE 
C 

WRITE(lZ>e000) TLABELClIf (ALA8ELm»I-l»3) 

6000 F0RHAT(63X*A10///7X*RE6RESSI0N>ESTINATE MSE DATA OVER 200 TRINOHIA 
IL SIMULATIONS. TWO REPLICATIONS PER CELL. DESIGN 2. P2A10*A2//) 

WRITE(12>8001I 

8001 F0RNAT(47X«A. ORIGINAL PRIOR IN BAYESIAN ESTIMATORS .«//// ) 
WRITEaZfBOlO) 

8010 F0RMAT(12X#27HPP****PP*P**P******P*P*P*P »*SS-25*»27H *♦*♦*♦♦♦♦♦** 
1****«4>*4<*«****,4X»28H ***********0************** »PSS-50P»27H 0000 
20000000000000000000000/0 ESTI- dir. 0,9H00*000000»0 PID-I5 *»10HP* 
3000^0t00,5Xf9H*********f* PI 0-40 ♦»10H*P**PPP**P*5X»9H*P**P****»* 
4PID-15 *#10HP*PPP*****»5X»9HP**P**P*P»* PID-40 ♦, 10HP****P****/* M 
5AT0R GEN. REPLICATION! REPLICATI0N2 REPLICATION! REPLICATI 
60N2 REPLICATION! REPLICATI0N2 REPLICATION! REPLICATION 

72*/) 

WRITE(!2»80!I) ( I» ( ( QLMSUf I» J>K )»K-!> 2)» J-!»2 ) > ( (QLMS12 ( I* J*K )>K 

!>!«2)*J-!»2)»I-!*!0) 

80!! FORMAT!* APMRO *I2» 4 (2E1 5.6« 2X ) / (9X»4f 2E!5.6>2X )) ) 
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WRITE 112*6012) ( I* < ( QLNS21 f I* J*K )»K>1*2 )* J*l* 2 )» ( ( QLNS22( I* J*K)* K 

1*1*2)*J-1*2)*I-1*10) 

6012 F0RHAT(« PHORO « 12* 4 ( 2E1 9.6* 2X I / < 9X* A ( 2E 15 .6*2 X I ) ) 

WRITE (12* 8013) (I* ( ( QLNS31 ( I* J »K ) *K«1 *2 ) * J-1* 2) * ( ( QLNS32( I* J*K) *K 

1-1*2)*J-1*2)*I-1*10) 

6013 FORNAT(* MLE * 12* 4 ( 2E15.6* 2X ) / (9X* A ( 2E15.6* 2X ) ) ) 

WR1TE(12*4997) 

WRITE(12*8000) TL ABEL ( 1 ) * ( AL ABEL ( I ) * I •! » 3 ) 

WRITE(12*8020) 

6020 F0RHAT(40X*B. UNIFORM AND PERTURBED PRIOR IN BAYESIAN ESTIMATORS. 
I*////) 

WRlTE(12*e010) 

WRITE (12*6021) (I* ( ( 0LMS41U* J»K) *K-1> 2) * J*l* 2)> ( ( QLMS42 ( I* J*K)* K 

1«1*2)*J*1.2)*1-1*10) 

8021 FORMAT(* APMRl P 1 2* 4 ( 2 E15.6* 2X ) / ( 9X* 4 ( 2E15 . 6* 2X ) ) ) 

WRITE (12* 8022) ( I* ( (QLMS5U I* J»K ) * K>1* 2 ) * J>1* 2 ) * ( ( QLMS52 (1 * J*K)*K 

1«1*2)*J-1*2)*I-1*10) 

6022 FORMAT(* APMR2 *I 2* 4 ( 2E1 5.6*2X) /( 9X*4 (2E15. 6»2X ) ) ) 

WRITE (12 *802 3) (I* ( ( 0LMS61 ( I* J *K ) *K>1* 2 ) * J«l* 2 ) * ( (0LMS62 ( I * J*K) »K 

1>1»2)*J*1»2)*I*1>10) 

8023 FORMAT)* PMDR2 *12* 4( 2E15.6* 2X) / t9X* 4 ( 2E15. 6*2X ) ) ) 

C 

9910 CONTINUE 

CALCULATE TUKEY DATA SUMMARIES* MEAN* AND STANDARD ERROR (S.E.) 

00 9925 IPI0«1»2 
DO 9925 NREPLIC-1*2 
SUM>0. 

SOSM-0. 

DO 9911 1«1*10 

SUM • SUN40LMS11(I*IPID*NREPLIC) 

SQSM-S0SM«-0LMS11( I* IPI0*NREPLIC)*0LMS11(I*IPI0*NREPLIC) 

9911 TUKEY(I)«0LMS11(I*IPID*NREPLIC) 

CALL SUMMARY) TUKEY* 10* Til (IPID*NREPLIC*1)*T11(IPID*NREPLIC* 2) *T11( 
1IPI0»NREPLIC*3)*T11(IPID*NREPLIC*4)*T11(IPID*NREPLIC*5) ) 
T11(IPI0*NREPLIC*6)-SUM/10. 

T11(IP10*NREPLIC*7)-SQRT((SOSM-10.*T11(IPID*MREPLIC*6)*T11(IPID*MR 

lEPLIC*6))/90.) 

T11(IPI0*NREPLIC*8)>T11( IPID*NREPLIC*3) 

T11(IPI0,NREPLIC*9)«(T11(IPID»MREPLIC*2)+2.*T11(IPI0*NREPLIC*3)4TI 
11(IPID*NREPLIC»4) )/4. 

T11(IPID*NREPLIC*10)-T11(IPID*MREPLIC»4)-T11(IP1D»NREPLIC*2) 

TliaPI0»NREPLIC*ll)«Tll(IPIO»NREPLIC*5)-Tll(IPID»NREPLIC»l) 

SUM«0. 
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SOSH-0. 

00 9912 

. SUM • SUH*0LMS12(I.IP1D»NREPLIC» 

S0Sn-SQSM«0LNS12(I*IPI0.NREPLIC}*0LMS12n*IPID»NREPLICI 

9912 TUKEYCI )-0LMS12U»IPI0»NREPLlCI 

CALL SUMMARY(TUKEY»10»T12{IP10»NREPLIC»1)»T12( IPID»NREPLIC»2I»T12I 
llPl0»NREPLlC»3»»T12nPn)»NREPLlC»<i)»T12tlP10»NREPLlC»5>) 
T12{IPI0»NREPLIC»6)»SUM/10. 

T12CIP1D»NREPLIC*7»-SQRT((SQSM-10.*T12IIPID»NREPLIC»6)PT12(IPID*NR 

lEPLIC*6n/90.) 

T12IIPI0»NREPLIC»8»-T12? IPI0#NREPLIC»3I 

T12CIPl0»MREPLIC»9»-(T12(IPI0.NRePLIC»2»^2.*tl2(lpr0»NREPLIC»3>4Tl 
12IIPI0*NREPLIC»4) )/4. 

T12<IPl0»NREPLIC,10)«T12(IPI0»MREPLIC»4»-T12nPI0,NREPLIC#2) 

T12nPlD»NREPLIC»ll)-Tl2(IPID»NREPLIC»5)-T12(rP10,NREPLJC*l) 

SUM-0. 

SQSM-0. 

00 9913 I-l»10 

SUM • SUM40LMS21(I#rPl0»NREPLICI 

SQSM-SQSM4QLMS21(I»IPl0»MREPLICI*QLMS2in»IPt0>NREPLICI 

9913 TUKEY(n-QLMS21(I»IPlO»NREPLIC) 

CALL SUMMARY(TUKEY,10#T21(IPIO,NREPLIC»1)»T21( IPI0»NREPLIC»2)»T21( 
IIPI0»NREPLIC»3>.T21( IPI0.NREPLIC»<»I.T21<IPID»NREPLIC,5n 
T2HIPI0»NREPLIC»6)-SUM/10. 

T21(IPl0»NREPLIC.7)-S0RT((SQSM-10.*T2inPlD.NREPLIC»6)*T21llPID»NR 

1EPLIC»6»)/90.) 

T2l(IPl0»NREPLrC»8l-T21(IPI0»NREPLIC»3) 

T21(IPl0f NREPLIC»9I-(T21{IPI0»NREPLIC»2|42.*T21(IPID»NREPLIC»3)4T2 
11(IPI0^NREPLIC.9) ) 

T21<IPl0»NREPLIC»10)-T21(IPID»NREPLIC»A)-T21(rPI0#NREPLIC»2) 

T2l<IPlO»NREPLIC»ll)-T2inPIO»NREPLIC»5)-T21(IPIO,NREPLIC»l> 

SUM-0. 

SQSM-0. 

00 9914 I-l»10 

SUM - SUM40LMS22(I. IPID.NREPLIC) 

SOSH-SOSM40LHS22(I»IPID.NREPLIC>*OLMS22(I»IPIO»MREPL1C> 

9914 TUKEY(n-QLMS22(I»IPI0. NREPLIC) 

CALL SUMMARY(TUKEY,10»T22(IPI0,NREPLIC.1)»T22(IPID,NREPLIC»2)»T22< 
llPl0»NREPLlC»3)»T22(IPI0»NREPLIC,4WT22nPI0»NREPLIC»5» ) 
T22CIPI0.NREPLIC»6I -SUM/10, 

T22(IPl0*NREPLlC»7)-SQRT((S0SM-lO.«T22aPID>NREPLIC*6}4T22(IP10»NR 

lEPLIC*6>)/90.) 

T22(IPI0>NREPLIC»8)-T22( IPI0>NREPLIC.3) 

T22(IPlD*NREPLIC>9)-(T22(IPI0»NREPLIC»2)42.«T22f IPI0>NREPLIC>3|4T2 
12(lPI0#NREPLIC»4))/4. 
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T22IIPI0»NREPLIC»10)-T22<IPID»NREPHC»<»)-T22(IPI0»NREPLIC»2) 

T22IIPI0»NREPLIC»11)-T22(IPI0»NREPLIC»5>-T22IIPI0»NREPLIC,1> 

sun-0. 

sosn>o. 

00 9915 I-l>10 ^ 

sun > SUn>QLnS31(I»IPI0»NREPLICI > 

SQSn-SOSn40LnS31(I*IPIO.NREPI.IC)*OLnS31(I>IPID*NREPLIC) 

9915 TUKEYCl )-0LHS31(I«lPID*NREPLIC) 

CALL SUMnARY(TUKEY»10»T31(IPI0.NREPLICf 1)*T31(IPI0#NREPLIC»2)«T31I 
1IP10«NREPLIC»3I*T31( IPI0*NREPLIC><t)»T31(IPID>NREPLIC»5) ) 

T31(IPI0.NREPLIC«6)- SUM/10. 

T31(IPIO»NREPLlC>7)-SORT((SOSn-10.*T31(lPID»NREPLlC>6)*T31llPlD»NR 

1EPLIC*6I)/90.I 

T31(1PID*NREPLIC«8)>T31(IPI0*NREPLIC*3> 

T31<1PI0»NREPLIC.9I« (T3UIPI0»NREPLIC»2)*2.PT31(IP1D,NREPL1C.3)*T3 
UUPlDfNREPLIC »<>))/ A. 

T31CIPI0»NREPLK»10)-T3UIPI0»NREPL1C»A)-T3HIPI0»NREPIIC»2» 

T3l<lP10*NREPLIC»lU-T31(IPI0»NREPLIC#5)-T31tIPI0»NREPLIC»l» 

sun-0* 

sosn-0. 

00 9916 l-lflo 

sun • SUn«OLMS32(I>IPID*NREPLlC) 

SOSn-SQSM*OLHS32(r» IPr0»NREPLIC)*QLHS32(I»IPI0.MREPLlC» 

9916 TUKEY(I)-QLMS3Z(I»IPI0,NREPLIC) 

CALL SUHMARY(TUKEY» 10» T3 2(1 PIO, MREPL IC» 1 1 » T32 ( I PI 0, MREPL IC» 2 » » T32 ( 
lIPIOfNREPLICf 3}>T32 ( IPI0*NREPLIC>A)»T32(IPID.NREPLIC»5) ) 
T32IIPI0,NREPLIC»6)-SUM/10. 

T32(IPI0«MREPLIC>7)-SQRT((SQSn-10.«T32(lPID«NREPLlC»6)«T32(lPID»NR 

XEPLlC»6l)/90.) 

T32(IPI0,NREPLIC»8)-T3Z( IPID#NREPLIC»3) 

T32(lPI0.NREPLIC>9)-(T32(IPI0.NREPLIC*2)«2.*T32(IPI0fMREPLIC»3)«T3 

12(IPID*NREPLIC»<in/A. 

T32(IPI0fNREPLIC»10l-T32(IPI0»NREPLIC*A)-T32(IPlD>NREPLIC>2) 

T32(IPI0»NREPLIC»11»-T32(IPI0»NREPLIC»5»-T32CIPID»NREPLIC»1) 

SUM- 0 . 

SQSn- 0 . 

00 9917 I-l»10 

SUM • SUM^QLMSAKI.lPIOfNREPLlC) 

S0SM-SQSM*0LMSA1(I,IPI0»NREPIIC1*0LHSA11I»IPID,NREPLIC) 

9917 TUKEY(I)-QLMS<*1(I»1PI0»NREPLIC) 

CALL SUMMARYdUKEYf 10fTAl(IPID.NREPLIC»l)>TAl(IPlO»NREPLlC>2)*TAl( 

. 1IPI0«NREPLIC»3)»T61(IPI0»NREPLIC»A}>TA1( IPID*NREPLICf5l} 
TA1(IPID»NREPLIC»6)-SUM/10. 

TAI(1PID»NREPLIC»71-SQRT((SQSM-10.*TA1(IPIO»NREPLIC*6»*TA1(IPID»NR 

lEPLIC>6))/90.) 



32 


T«KIPI0>NRE(>LIC>e)>T^l(IPI0*HREPLIC#3) 

T4l«IPI0»NREPLIC»9»»(T<*l(IRI0»NREPLIC»2l^2.*T^l(H>ID»NREPLIC»3»»T4 
llflPI0«NREPLIC>4) )/4. 

T41fIPI0«NREPLIC»10)«T41(IPID>NREPLIC*4l>T41(IPI0»NREPLIC»2) 

T^lfIPID»NREPLIC«ll)-T41(IPI0>NREPLIC#5)-T41(IPID>NREPLIC«l) 

SUN-0. 

SQSN-0. 

00 9916 I-l.lO 

SUN > SUN«QLNS^2(I*IPI0*NREPLIC) 

S0SN-SQSNi>0LNS42<I. IPI0>NREPLIC)«QLNS42(t*IPI0«NREPLIC) 

9918 TUKEY{II>QLNS42(I.IPI0«NREPLIC) 

CALL SUNNARY(TUKEY«10*T<)2(IPID»NREPLlC>l)»T^2tIPID»NREPLIC>2)*TA2( 
1IPI0«NREPLIC*3)»T42(IPI0»NREPLIC*4)«T42(IPI0.NREPLIC*5> ) 
T42(IP10*NREPLIC»6)-SUH/10. 

T42(IP10»MREPLIC»7)-SQRT((SQSM-10.*T%2CIPID»NREPLIC»6)*T42CIPI0*MR 

lEPLIC#6i)/90.) 

T^2<IPID#NREPLIC»8I-T42(IPI0,NREPLIC»3) 

T42(IPl0»NREPLIC»9I-<T<»2nPI0»NREPLIC#2>^2.*TA2nPID»NREPLIC»3> + TA 
12(IPI0»NREPLIC»4» 

T42(IPrO»NREPLICf 10l-T42(IPID#NREPLICf AI-T42(IPI0»NREPLIC»2) 
T42aPI0»MREPLIC»lll-T42(IPI0,NR£PLIC»51-T42(IPlO»NREPLIC»ll 
SUH-0. 

SOSN-0. 

00 9919 I-l.lO 

SUN • SUN^OLNSSKI.IPIO.NREPLIC) 

SQSN-S0SN+QLMS51(I»IPI0.NREPLIC»*QLNS51(I»IPI0#NREPLIC» 

9919 TUKEYCI »-0LNS51(I»IPID.NR£PLIC» 

CALL SUMMARY(TUKeY.lO.T5l(IPIO.NREPLIC»l).T51(IPIO»MREPLIC»21*T51( 
1JPID»NREPLIC»3I»TS1( IPI0.NREPLIC>4l»T51(IPID>NREPLlC»Sn 
T5iaPI0»NREPLIC»6)-SUM/10. 

T51(IPID»NREPLIC»7»-S0RT((SQSM-10.*T51(IPID»NREPLIC»6»*T51(IPID»NR 

1EPLIC.6II/90.I 

T51(IPID,NREPLIC.8I-T51(IPI0»NREPLIC»3) 

T5UIPID»NREPLIC»9)-(T51(IPI0»NREPLIC»2»+2.pT51UPID»NREPLIC#3)*T5 

ll(IPID*NREPLIC>4))/4. 

T51(IPID,NREPLIC.10l-T51(IPI0.MREPLIC.4»-T51tIPI0»NREPLIC»21 

T51(IPI0*NREPLIC»11)-T51(IPI0*NREPLIC»5)>T51CIPI0*NREPLIC»1} 

SUN-0. 

SQSN-0. 

DO 9920 I-l.lO 

SUN - SUN^QLNSS2(I»IPI0f NREPLIC) 

SQSN-SQSN«0LNS52(I«IPID»NREPLICI*QLNS52(I*IPI0>NREPLIC> 

9920 TUKEY(I)-0LHS52(I>IPID.NREPL1C> 

CALL SUNMARY(TUKEY»lO»T52<IPIO»NREPLXC»l)»T52(IPlO*NREPLXC»2l»T52C 
IIPI0*NREPLXC*3)>T52(IP10.NREPLIC>4)»T52(IPXD*NREPLIC*5) ) 
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T52(IPID»NREPLIC*6»«SUM/10. 

T52(IPI0*NREPLIC,7»«SQRT<(S0SM-lO.*T52(IPI0»MREPLIC»6)*T52(IPI0fNR 

1EPLIC*6>W90.I 

TS2nPI0»NREPLIC>8l«T52( IPI0«NREPLIC>3) 

T52CIPID»NREPLIC,«))«(T52(IPID»NREPLIC#2)*2.*T52HPI0»NREPLIC»3)+T5 
12(IPI0»NREPLIC>4) )/4. 

T52(IPI0*NREPLIC>10)-T52(IP10*NREPLIC*4)*T52(IPID»NREPLIC«2) 

T92(IPI0>NREPL1C«11)>T52(IPID«NREPLIC*5)-T52(IPID>NREPLIC>1) 

sun>o. 

SQSn-0. 

DO 9921 

SUN • SUM«QLMS61(I*IPI0»NREPLIC} 

S0SH-S0SN«0LNS61(I>IPI0*NREPLICI*0LNS61(I»IPID»NREPLIC) 

9921 TUKEY<n- 0 LHS 61 (I»IPID.NREPLICI 

CALI SUMMARY(TUKEY#10»T61(IPI0#NREPLIC#1)#T61(IPI0>NREPLIC»2»»T61C 
1IPI0>NREPLIC#3)>T61< IPI0«NREPLIC>AI»T61(IPID>NREPLIC*9M 
T6l(IPI0»NREPLIC»6)«SUM/10. 

T61<IPI0»NREPLIC»7l-SQRT<(SQSN-lO.*T61(IPID*NREPLIC>6)9T61(IPID#NR 

lEPLIC»6l)/90.) 

T61(IP10»NREPLIC*ei>T61(IPIO»NREPlIC»3) 

T61(IPI0#NREPLIC»9)-{T6inPID»NREPLlC»2)+2.*T61( IPID»NREPLIC»3»*T6 
linPIO*NREPLIC»An/A. 

T61(IPI0»NREPLIC»10»-T61(IPI0»NREPLIC»^)-T61UPIDf NREPIIC»2» 
T61(IPI0»NREPLIC»ll)-T61(IPI0»NREPLIC»5l-T61(lPID»NREPLIC»l) 

SUN-0. 

SOSN-0. 

DO 9922 I-lflO 

SUN > SUN^0LNS62(I.IP10»NREPLIC) 

S0SN-S0SN*QLNS62(I«IPI0»NREPLICI*QLNS62(I»IPID»NREPLIC) 

9922 TUKEY<n-QLHS62(I»IPI0»NREPLIC> 

CALL SUNNARY(TUKEY»10»T62(IPID>NREPLIC»1)>T62< IPID«NREPLIC»2I>T62( 
llPIOf NREPLIC»3)»T62(IPID.NREPLIC»9}.T62f IPID»NREPLIC»5) ) 
T62(IPI0»NREPLIC»6»"SUM/I0. 

T62(IPI0>NREPLIC>7)-S0RT<(SQSN-lO.«T62(IPI0»NREPLIC*6)*T62llPID*NR 

lEPLIC»6l)/90.) 

T62(IPI0>NREPLIC>8I-T62(IPI0»NREPLIC>3) 

T62(IPI0>NREPL1C.9)>(T62(IP10»NREPLIC>2)^2.AT62(IPID»NREPL1C»3)*T6 
12(IPI0»NREPLIC»A) lYA. 

T62IIPI0»NREPLIC»10)«T62(IPI0»NREPLIC»AI-T62CIPI0»NREPLIC»2» 
T62(IPI0>NREPLIC>11)-T62f IPI0»NREPL1C>$)-T62(IPI0*NREPLIC«1) 

9925 CONTINUE 
NREPLIC-1 

UR1TE(9»5100) (ALABELn)*I*l»3) 

5100 F0RNAT(63X«TABLE 7.5«///* DATA SUNNARIES* CENTRAL VALUES* AND SPRE 
IADS* NULT. BY 10* OVER 10 OIRICHLET SINULATIONS FOR Q.L. REG-EST N 



2SE. FIRST REPLIC* DESIGN 2.«///57X*3A10/////> 

4990 NSS-25 
PID-0.15 
ITABLE-0 
WRITE(5»5105) 

5105 F0RHATa9Xt53H****P******P********DATA SUHMAR Y******************* 
1A5X»36H********P**CENTRAL VALUES*********P*4X17H*****S PR EADS *♦*•♦/ 
2* ESTIMATOR SS PIO L EXTREME L HINGE MEDIAN U HINGE U 
3EXTREME*7X«MEAN (S.E.) MEDIAN TRIMEAN MIDSPREAD RANGE* 
4/1 

4991 WRITE (5* 5510 1 NSS* PIO* ( TIK IPIO*NREPLIC»K 11 ) > (T21 ( IPIDvNREP 
UIC*K)*K>1*11)* (T31(IPI0*NREPLIC*K)>K>1*11)» (T4K IPIO*NREPLIC»K)*K 
2>l*lll* (TS1(IPI0*NREPLIC»K)*K-1*11)* ( T61 ( I PI 0*NREPL IC * K ) > K>1* 11) 

5510 F0RMAT(3X*APMR0*3X»I2»F4.2»5F11.7*3XF9.6*«*F9.7*)*2F9.6»3X*2F9.6/3 
1X*PMDR0*9X*5F11.7*3X.F9.6*(*F9.7*)*2F9.6»3X»2F9,6/3X*MLE*11X*5F11. 
27»3X»F9«6*(*F9.7*)*2F9,6»3X*2F9.6//3X*APMR1*9X»5F11,7»3X*F9,6*(*F9 
3.7*I*2F9.6*3X*2F9.6//3X*APMR2*9X»5F11.7»3X»F9.6*<*F9.7*)*2F9.6*3X» 
42F9,6/3X*PM0R2*9X»5F11.7»3X»F9.6*I*F9.7*)*2F9.6,3X*2F9.6/) 
ITABLE-ITABIE*! 

60 TO <4992*4993*4994*4995) ITABLE 

4992 PI0«0«40 
60 TO 4991 

4993 NSS-50 
PIO-0.15 
60 TO 4991 

4994 PID-0.40 
60 TO 4991 

4995 WRITE(5*5112) 

5112 F0RMAT(///3X* NOTE THAT A ZERO BEFORE A DECIMAL DENOTES AN EXACT Z 
lERO. OTHERWISE* THE ZERO IS ROUNDED.*) 

IF (NREPLIC-1) 4996*4996*4998 

4996 NREPLIC-2 
WRITE(5*4997) 

4997 FORMAT!/////) 

WRITE(5*5115) <ALABEL<I)*I-1*3) 

5115 FORMAT(63X*TA6LE 7.6*///* DATA SUMMARIES* CENTRAL VALUES* AND SPRE 
IADS* MULT. BY 10* OVER 10 DIRICHLET SIMULATIONS FOR Q.L. REG>EST M 
2SE. SECOND REPLIC. DESIGN 2.*///57X*3A10/////) 

60 TO 4990 

499B CONTINUE 
END 
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FUNCTION GAHHAtGG) 

GENERATE A GANNA RANOON VARIABLE 

TO 00 SO USE ALGQRITHN GT FRON AHRENS AND DIETER (197A) "CONPUTER 
NETHOOS FOR SANPLING FRON GANNA» BETA> POISSON* AND BINONIAL 
DISTRIBUTIONS"* P229* CONPUTING. VOL. 12 

NOTE THAT FOR 1/77 SINULATION STUDY* GG RANGES FRON 0.1 TO 9.8 

OBTAIN INTEGER PART OF 66 

K>6G 

OBTAIN FRACTIONAL PART OF GG 
F-GG-K 

OBTAIN INTEGER PART OF GANNA 

6I«0. 

IF (K-OI 19*14*8 
8 I-O 
GP"1. 

10 I>I^1 

UU«URAN(0.) 

GP-GPPUU 

IF (I>K1 10*12*12 
12 6I—AL0G(GPI 

OBTAIN FRACTIONAL PART OF GANNA 
19 6F»0. 

IF IF-0.) 90*90,15 

15 B«(2.7182818289590*F)/2.7182818269590 
OF-l./F 

FNINl-F-1. 

16 UU-URANIO.) 

6P«B*UU 

GENERATE NEW UNIFORN RANOON NUNBER FOR TESTS IN FOLLOWING STEPS 
18 AND 30 

UU>URAN(O.I 
IF (6P-1.I 18*18*30 



o o o 


C 6F IS LESS THAN OR EQUAL TO 1. 
C 

18 6F«GP**0F 
TEST-EXPC-GFI 
IF (UU>TESTI 40*40*16 

6F IS GREATER THAN 1. 

30 GF»AL0G((B-6P)*0F) 
TEST"GF**FMIN1 
IF (UU-TEST) 40*40*16 
40 GANNA-6I4GF 
RETURN 
ENO 
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C 

C 

SUBROUTINE 6ENXZ(UU>NSSI 


GENERATE NULTINQNIAL CONPLETE (X) AND INCOMPLETE IZ) DATA 
ALSO CALCULATE HAXINUN-L IKEL IHOOD-ESTIMATES FROM COMPLETE DATA 


COMMON OEP(3*3)«DQL(^*3)«E(4»3)»IROBUST>NTS>Pl»P2*P3»TPID 
COMMON/CALEST/APHClI*APHC12>APHC22»CONVCRI»COVSKIP»DMLCl»OMLCZ>DML 
1C3«OPIO*EAPMC11«EAPMC12»EAPMC22>EPMC11*EPHC12»EPMC22»ISTOP*N12»N13 
2»N23«PIO»PMLCl*PMLC2*PMLC3>SS#SSN»TlMAP#riMEP»riM(2)«TIN21*TIM31>X 
3NUl«XNU2«XNU3*Zl*Z2f Z3*Z12>Z13>Z23»Z1N>Z2N»Z3N 
C0NM0N/0ATA/XDATA(2)«Z0ATA(61 
DIMENSION UU(NSS) 

INTEGER Y1.Y2.W1*W3»V2#V3»C0VSKIP 

EOUIVALENCE (E t 1» 1 ) » PE PMl) . I E U» 2 » *PEPM2 > , ( E t 1» 3 » ♦ PEPM3 i 
(E(2>l)f PMLE1)>(E(2»2)»PMLE2)« (E(2»3WPMLE3) 
(E(3*1)*PPMD1)>(E<3*2I»PPMD2>> (E(3>3)*PPMD3) 
(E(A»1I»PAPM1)>(E(A»2)»PAPM2)* (E(4*3>»PAPH3) 
(OGP(l>li*EHLl)> (0EP(1*2>*EML2I> ( OEP ( 1* 3 )» EML3I 

(0EP(2«1)»EPM01)> f 0EP(2»2)*EPM02)» ( OEP (2» 3 ) > EPMD3) 
(OEP(3>l)tEAPMNl)» (0EP(3»2)»EAPMN2)>(0EP(3*3)>EAPMN3) 
(OOL( l>ll>DEPMNU*(D0L(l*Z)»DEPMN2)>(00La>3)»DEPMN3) 
(D0L(2»1)>0ML1)» (0QL(2*2)»DML2)f ( DDL ( 2» 3 ) » DML 3 ) 

(DaL(3»l)»0PMDl)» (DOL(3»2)»DPMD2)» (DQL ( 3> 3 ) » 0PM03) 
(DQL(A»I)»DAPNN1I» t0QL(^>2)»0APMN2)>(0QL(A>3)»0APMN3) 


EOUIVALENCE 

EOUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EOUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EOUIVALENCE 


RECALL THAT NSS IS THE INTEGER FORM OF THE SAMPLE SIZE SS 
CALCULATE APPROX AMOUNT OF DATA GOING INTO EACH OF THE 4 GROUPS 


HA-PID/2. 

H3-HA 

H2-H4 

H1»1.-3.*H2 

E3-H1+H2 


SET END POINTS 

P12«P1*P2 

E1-P1*E3 

E2-P12«E3 

EA>E3>P1*H2 

£5-E3*P12*H2 

E6*H1^Z.«HZ 
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E7«E6*P1*H2 

E8-E6«P12«H2 

INITIALIZE Z AND OUHNY VARIABLES Yl» UI> AND VI 

Z1>0. 

Z2-0. 

Z3-0. 

Yl-0 

Y2»0 

Wl-0 

W3-0 

V2«0 

V3>0 

GENERATE X* Z DATA 


CALL UNIFORN RSEUOO RANOQM«NUNBER GENERATOR 
XIN+11 • (A3A902756A7AA5.*X(Nn HOO( 2EXP J ) 

SPECTRAL NUMBERS C(2)*2.839# C(3)«2«095> C(AI>1.B19» Cf$)«0.97B 

CALL URANVIO.»NSS»UUI 
C 

00 85 I-1»NSS 
U«UU(I) 

IF (El-U) 2»2»A0 
2 IF (E2-UI A»A»A5 
4 IF IE3-U) 6»6»50 
6 IF (E4-U) 8>8*55 
8 IF (E5-U) 10»10»60 
10 IF (E6-U) 12»12»65 
12 IF CE7-U) 14»14»70 
14 IF IE8-U) 80»80»75 
40 zi-zm. 

GO TO 85 
45 22»22*1. 

GO TO 85 
30 Z3«Z3-»1. 

GO TO 85 
55 Yl-Yl+1 
GO TO 85 
60 Y2-Y2+1 
60 TO 85 
65 W3>W3«1 
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70 

GO TO 85 
75 V2-V2^1 
60 TO 85 
80 V3«V3«1 
85 CONTINUE 
N12-YI+Y2 
N13«W1>W3 
N23»V2fV3 

OBTAIN REAL FORM OF SHARED INCONPLETE DATA 

Z12«N12 

Z13-N13 

Z23-N23 

OBTAIN CQNRLETE DATA X 

Xl«Zl*YlfHl 

X2«Z24Y2^V2 

K3-Z3+W3+V3 

CALCULATE COMPLETE-DATA NAXIHUHH IKEL IHOOO ESTIMATES 

PMLCl-Xl/SS 
PHLC2-X2/SS 
PHLC3»1.-PMLC1-PHLC2 
IF (PMIC3-0.J 90f95»95 
90 ISTOP-1 

PRINT 92» XNUlf XNU2»XNU3»PI»P2»P3»PI0»SS»X1»X2»X3»Z1»Z2#Z3»Z12»Z13 
l*Z23»PHLClf PMLC2.PMLC3»NSS»NTS 

92 FORHATi///* GENXZ. COMPLETE-DATA HLE NEGATIVE P. XNU»*3F10,4»* GE 
INERATEO P-*3F10.6»* PID«*F6.2»* SS-*FA,0/* X-*3F6.0»* Z-*6F6.0»* P 
2HL«*3F10.6.* NSS-*I3»AX»I3* TH TRINOMIAL SINUL*) 

RETURN 

95 OMLCl-PMLCl-Pl 
0MLC2-PMIC2-P2 
0HLC3-PMLC3-P3 
ZlN-Zl+XNUl 
Z2N-Z2+XNU2 
Z3N>Z3«XNU3 

PUT X AND Z DATA INTO SEPARATE BLOCK COMMON BECAUSE PROGRAM WQN«T 
CORRECTLY RUN IF BLANK COMMON X>Z DATA IS PUT INTO KTITER* COUNTS* 
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C AND BESTEST (PERHAPS PROBLEMS WITH E EQUIVALENCE STATEMENTS) 
C 

XDATA(1)«X1 
X0ATA(2)-X2 
Z0ATA(1)-Z1 
ZOATACZI-ZZ 
ZDATA(3)<iZ3 
ZDATA(4)>Z12 
Z0ATA(5)-Z13 
Z0ATA(6)*Z23 

TRUE PERCENT (PROPORTION) INCOMPLETE DATA 

TPI0-(Z12«Z13«Z23)/SS 
OPIO-PIO-TPIO 
RETURN 
END 
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SUBROUTINE EPN 

RROGRAN TO CALCULATE EXACT POSTERIOR MEAN AND COVARIANCE MATRICES 

CQHNON aEP(3«3)»OOL(4»3}»E(4#3)f IR0BUST»NTS»P1*P2»P3*TPID 
COMMON/C ALEST/APMCH»APHC 12* APHC22»CONVCRI,COVSK I P»0HLC1»0NLC2»0NL 
IC3 pOPIO.EAPMC11pEAPMC12»EAPMC22»EPHC11»EPHC12»EPNC22»ISTOP»N12,N13 
2pN23pP10pPHLC1pPMLC2pPMLC3pSSpSSNpTIMAP,TIMEPpTIH(2)pTIM21pTIM31pX 
3NUlpXNU2pXNU3pZlpZ2pZ3pZ12pZ13pZ23priNpZ2NpZ3N 
INTEGER CQVSKIP 

EQUIVALENCE ( E ( Ip 1 1 p PEPNDp f E f 1 p2 )p PEPH 2 )p ( E( Ip 3) p PEPN3 ) 
EQUIVALENCE ( E( 2p 1) p PHLE 1 Ip ( E I 2p 2 1 p PNL E2 I p ( E (2p 3 ) p PNL E3 I 
EQUIVALENCE ( E< 3p 1 1 p PPHOl )p ( E< 3 p 2) p PPN02 Ip IE ( 3p 3 1 p PPND 3 I 
EQUIVALENCE ( E ( 4p 11 p PAPN l Ip (E( 4p 2 Ip PAPN 2 Ip I E ( Ap 3 I p P APN3 1 
EQUIVALENCE (OEP (Ip 1 1 p EHLl I p ( DEP ( Ip 2 I pENL2 I p ( DE P ( Ip 31 p ENL3 1 
EQUIVALENCE ( OEP ( 2p 1 1 p EPHOI I p ( OEP ( 2 p 2 I p EPHD 2 1 p ( OE P ( 2p 3 1 p E PN03 1 
EQUIVALENCE ( OEP ( 3p 1 1 p E APHNl |p ( DEP(3 p2 ) p E APHN2 ) p ( DE P ( 3p 3 1 p E APNN3I 
EQUIVALENCE ( OQL ( Ip 1 1 p OEPNNl I p ( DQL ( I p 2 I p 0EPNN2 Ip ( DQL ( Ip 3 1 p DEPNN 3I 
EQUIVALENCE (OOL (2p 1 1 p ONLI I p ( OOL (2 p 2 I p 0ML2 Ip ( OQL ( 2p 3 I p OML 3 I 
EQUIVALENCE ( OOL ( 3p I Ip OPNDI Ip ( OOL (3p 2 I p 0PM02 1 p (OOL ( 3p 3 I p 0PHD3 I 
EQUIVALENCE (OQL (^p 1 1 pOAPHNI |p ( OOL (<>p 2 I pDAPNN2 Ip ( OQL ( ^p 3 |p 0APNN3 1 

N121«N12>1 

N131«N13>1 

N231»N23>1 

3 Yl-ZIN 
YY1A»GAM(Y1I 

IF (ZIZ'O.I <>p<»p8 

4 IF (Z13-O.I 5p5p8 

5 IF (Z23-O.I 6p6p8 

CONPLETE-OATA CASEp ALL ZIJ>0. 

6 SUNI2>YY1A*GAN(Z2N|PGAN(Z3N| 

S12P1»Z1N*SUN12 

S12P1SQ-(Z1N+1,I*S12P1 

. S12P2»Z2N*SUHI2 
. S12P2SQ«(Z2N+1,I*S12P2 
S12P1P2-Z1N*Z2N*SUH12 
60 TO 55 
8 Z2N12-Z2N+Z12 

Z3N1323«Z3N+Z13*Z23 

SIJPK DENOTES ZIJ SUN FOR POSTERIOR MEAN P(K) CALCULATIONS. 
SINlLARLYp TIJPK DENOTES A TERN OF THIS SUN. 
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S12P1-0. 

S12P2-0. 

S12P1SQ-0. 

S12P2SQ-0. 

S12P1P2-0. 

SUN12-0. 

C 

00 50 1IA>1»N121 
IA-IIA«1 

IF ClA-Ol 9,q»10 
9 C0N12-1. 

Y2-Z2N12 

YY2A-6AH(Y2) 

GAtinY3-GAM<Z3N1323) 

GO TO 16 

10 CON12«(Z12-IA+1.)*CON12/IA 
YY1A-Y1*YY1A 
Y1«Z1N>IA 
Y2-Z2N12-IA 
YY2A-YY2A/Y2 

16 SUfU3«0. 

S13P1-0. 

S13P2-0* 

S13P1SQ-0. 

S13P2SQ-0. 

S13P1P2-0. 

ZINIA-ZIN^IA 

C 

00 60 IIB«1>N131 

IF (IB'OI 17«17>20 

17 C0N13-1. 

YYIB-YYIA 
Y3-Z3N1323 
YY3B-GAMHY3 
GO TO 27 

20 BI-IB-1 

CON13-(Z13-BI)«CON13/IB 
Y3-Z3N1323-IB 
YY3B-YY3B/Y3 
YY1B-YY1B*(YI>BII 
27 YY2C-YY2A 
. SUM23-0. 

S23P2-0. 



S23P2SQ-0 


00 39 IIC«1»N231 
IC-IlC-1 

IF (IC>0) 30»30*31 

30 C0N23-1. 

YY3C-YY3B 
00 TO 34 

31 Cl-IC-1 

CON23>(Z23-CI)«CON23/IC 

YY2C-<Y2^C1)*YY2C 

YY3C-YY3C/IY3-IC) 

34 T23«C0N23*YY2C*YY3C 
F«Y2*IC 
T23P2-T23*F 
T23P2SQ-T23P2*(F*1. ) 
SUH23-SUM23»T23 
S23P2-S23P2+T23P2 
S23P2SQ«S23P2SQ^T23P2SQ 

39 CONTINUE 

G«C0N13*YY18 
ee«ZlNU«IB 
T13"G*SUM23 
. T13P1«T13*GG 
T13P2»G*S23P2 
T13P1SQ"T13PI*<G6^1.I 
T13P2S0»G*S23P2S0 
T13P1P2»T13P2*GG 
SUH13»SUN13*T13 
.. S13P1-S13PI+T13PI 

S13P2-S13P2+T13P2 
Sl3PlS«l-S13PlSO+ri3PlSQ 
S13P2S0-S13P2S0*T13P2S0 
S13P1P2-S13P1P2*T13PIP2 

40 CONTINUE 

T12-C0N12*SUN13 
T12P1»C0N12*S13P1 
T12P2-C0N12*S13P2 
T12P1S0»CQN12*S13P1SQ 
. T12P2S0«C0N12*S13P2S0 

T12P1P2-C0N12*S13P1P2 
SUH12«T12>SUN12 
.. S12PI-S12P1 + T12P1 
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S12P2-S12P2>ri2P2 
S12P1SQ»S12PISQ>T12P1SQ 
■ S12P2S0»S12P2S0^T12P2SQ 
SI2P1P2-S12P1P2+T12P1P2 
50 CONTINUE 

ELEMENTS OF POSTERIOR MEAN OF P GIVEN Z 

55 01»SUM12*SSN 
02«0I*CSSN>1.) 

PEPM1-S12P1/01 
PEPM2-S12P2/01 
PEPM3-1.-PEPM1-PEPM2 
IF <PEPM3-0.) 230»201.20l 
201 CALL SEC0N0(TIM2I 

ELEMENTS OF POSTERIOR COVARIANCE MATRIX OF P GIVEN Z. 

XIX2-S12PIP2/02 

EPMC12"XIX2-PEPMI*PEPN2 

X1SQ«S12P1SO/02 

EPNC11«X1S0-PEPMI*PEPMI 

X2SQ«S12P2S0/D2 

EPMC22»X2SO-PEPH2*PEPN2 

OEPMNl-PEPHl-Pl 

0EPMN2-PEPM2-P2 

0EPMN3-PEPH3-P3 

RETURN 

230 PRINT 231 

231 FORMAT!////* EXACT POSTERIOR MEAN, NEGATIVE P.*//> 

PRINT 2A0» PEPMl,EPNCll»XlSQ»Sl2Pl»S12PiSQ»PEPM2»EPMC22*X2SQ*S12P2 
1>S12P2SQ»PEPM3»EPHC12»X1X2»SUN12»SI2P1P2*NTS 
2A0 FORMAT! //7H PEPMl«El5,8»AXf 8H EPHCll«El5.8f AX*6H X1SQ-E15.8* AX*7H 
1S12P1«E15.8»AX»9H S12P1S0-E15.8/7H PEPM2-E15.8* AX»8H EPMC22-E15.8, 
2AX»6H X2SQ>E15,8»AX»7H SI2P2-E15.8>AX,9H S12P2SQ-E1S.8/7H PEPH3«EI 
35,8«AX»8H EPMC12>E15.8fAX»6H X1X2«E1A.7» AX»7H SUM12*E1A,7*AX»9H SI 
A2P1P2-E1A.7** NTS-*I3/» 
lSTOP-1 
RETURN 
END 
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ORIGINAL I?A(JE 
OF POOR QUALOT 


FUNCTION GANtXI 

EVALUATE GANfXI FOR CASES FROM PRIOR C.I».l>9.8» 
IF (X-9.21 1,35*35 


ARGUMENT IS LESS THAN OR EQUAL TO 9.1. TO INSURE 11 SIGNIFICANT 
PLACES OF ACCURACY IN GANNA, USE 11 SIGNIFICANT-FIGURE VALUE 
CALCULATED FROM LCG(GANIXI) FROM ABRANOUITZ AND STEGUN OR DAVIS. 
CALCULATE GANIO.ll, GANll.l), AND 6ANI2.1) BY HAND FROM 6ANI3.1) 
AND RELATION GAN(X«1I> X GAHIX} 


C 


C 


C 


c 


c 


c 


c 


c 


1 IF (ABS(X-9.11-1.E-13I 2,2, A 
GANI9.il 

2 6AN-A9973.70B949629 
RETURN 

A IF (ABS(X-B.1)-1.E-13) 6,6,8 
GANIB.ll 

6 GAN-6169. 593697A851 
RETURN 

6 IF (ABS(X-7.1I-1.E-13I 10*10,12 
GANI7.il 

10 GAN-868.95685680072 
RETURN 

12 IF IABS(X-6.1I-1.E-13I 1A,1A,16 
GAN(6.ll 

lA 6AN-IA2.A519AA06569 
RETURN 

16 IF (ABS(X-5.1)-1.E-13I 18,18,20 
GANC5.il 

18 GAN-27.931753738371 
RETURN 

20 IF (ABS(X-A.1I-1.E-13I 22,22, 2A 
6ANCA.il 

22 GAN-6. 8126228630175 
R E TUR N 

2A IF CABSCX-3.1I-1.E-13I 26,26,28 
6ANC3.il 

26 GAN-2.1976202783927 
RETURN 

28 IF CABSCX-2.1I-1.E'13I 30,30,32 
6ANC2.il 

30 GAN-1. 0A6A858A685A 
RETURN 
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32 IF (ABS(X>l.ll-l.E>13t 34»3A»36 
C 6AN(1.1I 

3A GAn>0.9513S076987 
RETURN 

C GANCO.il 

36 GAN-9.513S076987 
RETURN 

ARGUNENT X IS GREATER THAN OR EQUAL TO 10. 

USE STIRLING«S FORMULA TO OBTAIN AN APPROXIMATION TO GANNA 
THAT IS ACCURATE TO 11 SIGNIFICANT FIGURES 

35 XSQ-XPX 
XCU-XSOPX 
XFIFTH-XSOPXCU 

Y IS THE APPRQXINATEO NATURAL LOGIBASE El OF GANNAIXI 

Y«(X-0.5l*ALOG<XI-X*0.91893853320A67*l./(12.PXI-l./C360.*XCUI*l./C 
11260. ♦XFIFTHI 
IF CX-22.1 A0»A5»A5 
40 Y«Y-1./(1680.PXFIFTH*XS0I 
45 GAN-EXPCYI 
RETURN 
ENO 


\ 
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SUBROUTINE NETHQOS 

FOR INCONRLETE DATA CALCULATE NY APPROXIMATION* POSTERIOR MODE* 

AND MAXIMUN-LIKELIHOOO ESTIMATE 

INTEGER COVSKIP 

DIMENSION A(3>3)«B(3»l}»IPIVOTt3)*lWKf6) 

COMMON OEP(3>3)*OQL(4>3t*E(A*3>*IROBUST»NTS*Pl»P2*P3*TPID 
COMMON/C ALEST/APHC 11* APMC12*APMC22*C0NVCR I »CaVSKIP*0MLCl*0MLC2*0ML 
IC3»OPIO»EAPMC11»EAPMC12»EAPHC22»EPMC11,EPMC12.EPMC22»ISTOP»N12*N13 
2»N23*PID«PMLC1*PMLC2*PMLC3»SS*SSN*TIMAP*TIMEP*TIM(2)*TIM21*TIM31*X 
3NU1*XNU2»XNU3*Z1>Z2*Z3*Z12*Z13»Z23>Z1N»Z2N*Z3N 
COMMON/ ITKT/AVNUMIT (6 »*CTMUMIT<6»7> 

EQUIVALENCE ( E ( 1» 1) * PEPMlt* ( EC 1* 2 I > PEPM2 ) * (E (1*3 ) * PEPM3) 
EQUIVALENCE ( E < 2* 1 ) • PMLEU* ( E (2*2 ) * PMLE2 ) * ( E (2*3) * PHLE3 ) 
EQUIVALENCE (E (3* 1 ) * PPMOl )* ( E(3* 2 )* PPMD2 ) » (E (3*3 > * PPH03) 
EQUIVALENCE (E ( A* I) » PAPMD* (E ( A»2 )• PAPM2 ) * ( E( A* 3)* P APM3 ) 
EQUIVALENCE (OEP(l*l)*EMLD* ( OEP( 1*2 >* ENL2 )* (OEP ( 1* 3 )*EML3) 
EQUIVALENCE (OEP (2* 1)»EPN01)» ( OEP ( 2*2 ) * EPMD2 > * (0EP(2»3I*EPM03) 
EQUIVALENCE (OE P ( 3* II » EAPMND* ( OEP 0*2 ) * EAPMN2 I* (OEP ( 3* 3) * EAPMN3) 
EQUIVALENCE (OOL ( 1* U * OEPMND* ( OQL ( 1*2 1* 0EPNN2 )* ( OQL ( 1* 3 ) *0EPMN3 ) 
EQUIVALENCE (0QL(2*1 )*OMLD* ( 0QL(2*2>* DML2I* (0QL(2*3)*DML3) 
EQUIVALENCE ( OOL ( 3* 1 I *0PM01I* ( OOL ( 3*2 ) * 0PM02 I* (OQL ( 3* 3 ) *DPM03 I 
EQUIVALENCE (OQL( A* 1) *0 APMNII * ( OQL (^*2 1 * 0APMN2 1* (DQL(4*3I»DAPMN3) 

CHECK TO INSURE THAT CONVERGENCE CRITERION IS NOT TOO STRICT 

IF (IROBUST-1) 1*100*90 

INCOMPLETE-OATA MAXIMUM-LIKELIHOOD ESTIMATE 

1 PMLEl-PEPMl 
PMLE2-PEPM2 
PMLE3-PEPM3 
K-1 . 

5 PLl-PMLEl 
PL2-PMLE2 

PMLE12-PMLE1+PMLE2 
PMLE13-PMLEI+PMLE3 
PMLE23-PMLE2+PMLE3 
IF (PMLE12-l.E-l<») 7,7*8 

7 TEMP-0. 

CO TO 9 

8 TBMP-Z12/PMLE12 

9 PMLEl-(Zl*PMLElP(TEMP*Zl3/PMLE13n/SS 
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PHLE 2 -lZ 2 ^»>MLE 2 *CTEHP 4 Z 23 /l»HlE 23 n/SS 

PHLE3-1.-PMLE1-PHLE2 

IF (PMLE3>0o) 21»l^tl4 

14 IF CPMLEl-O.n 19*15*16 

15 IF (ABS(PHLE1-PL1I*0. 00001) 17*20*20 

16 IF (ABS(PMLEl-PLl)/PflLEl-CONVCftl) 17*20*20 

17 IF IPMLE2-0.1) 18*18.19 

18 IF (ABS(PMLE2-PL2)<-0. 00001) 25*20*20 

19 IF (ABS(PNLE2-PL2)/Pf1LE2>CQNVCRI) 25*20*20 

20 K-K^l 

IF (K'lOOO) 5*5*21 

21 PRINT 22* K*XNU1*P1*PL1*PL2*C0NVCRI*IR0BUST*NTS*TPID*PNLE3»XNU2*P2 
1*PNLE1*PHLE2 

22 FORNATC// * EXCESS NUMBER OF ITERATIONS FOR INCOMPLETE-DATA M.L.E 
1* 1S*I3»* XNU1**F9,4.* P1-*F9.6*3X*PL1-<!F11,8*3XPPL2»*F11.8/* COMV 
2CR0)«*F7,5* IR0BUST-PI2* NTS»*I3* TPID«*F6.2* PMLE3-PF6.4* XNU2»*F 
39.49 P2-9F9.6* PMLE 1«*F11.8* PMLE2-*F11.8 ) 

IF (K-1030) 25*25*250 

CONVERGENCE FOR MAXIMUM-LIKELIHOOO ESTIMATE INCOMPLETE DATA' 

2S EML1«PMLE1-PEPM1 
EML2-PMLE2-PEPM2 
EHL3-PMLE3-PEPM3 
0HLI"PMLE1-P1 
0ML2-PMLE2-P2 
> DH13-PMLE3-P3 

INCREMENT COUNTER FOR NUMBER OF ITERATIONS FOR CONVERGENCE 
CALL KTITERIK.I) 


POSTERIOR MODE 

SO Tl-ZlN-1. 
T2-Z2N-1. 

K-1 

. O-SSN-3. 
PPMOl-PEPMl 
PPM02-PEPM2 
PPMD3-PEPM3 
.55 PLl-PPMOl 
>. PL2-PPMD2 
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PPN012-PPHD1+PPM02 
PPH013-PPH01 + PPM03 
PPM023-PPM024PPM03 
IF IPPM012-1.E-14) 57*57*58 

57 TEMP-0. 

60 TO 59 

58 TEMP-Z12/PPMD12 

59 PPM01-(Tl*PPMDl*<TEMP^Z13/PPM013n/0 
IF (PPM01 .lt. 0. ) PPMOl-0. 

PPM02-(T2*PPM02*(TEMP*Z23/PPM02 3n/0 
IF (PPM02.LT.0. ) PPM02-0. 

PPM03-1.-PPM01-PPM02 

IF (PPN03-0.) 71*64*64 

64 IF (PPMOl-0.1) 65*65*66 

65 IF (ABS(PPMOl-PLl)-O.OOOOll 67*70*70 

66 IF (ABSIPPMDl-PLll/PPMDl-CONVCPn 67*70*70 

67 IF (PPM02-0.il 68*68,69 

68 IF (ABS(PPMD2-PL2)-0. 00001) 75,70*70 

69 IF (ABS(PPM02-PL2)/PPMD2-C0NVCRI) 75*70*70 

70 K-K+l 

IF (K-lOOOl 55,55*71 

71 PRINT 72* K*XNU1*P1*PL1*PL2*CONVCRI*ZROBUST*NTS*TPIO*PPM03»XNU2*P2 
1*PPMD1*PPH02 

72 FORMAT!// ♦ EXCESS NUMBER OF ITERATIONS FOR INCOMPLETE-DATA PPMO 
1. IS*I3,4 XNU1-*F9.4** Pl-*F9.6*3X*PLl-*F11.8*3X*PL2-*Fll.e/4 CONV 
2(RD)-*F7.5* IR0BUST-*I2* NTS-*I3^ TPID-*F6.2* PPH03-*F6.4* XNU2-*F 

-39. 4* P2-*F9.6* PPM01-*F11.8* PPM02-*F11 . 8 ) 

IF (K-10301 75*75*250 

CONVERGENCE FOR POSTERIOR MODE INCOMPLETE DATA 

75 EPMDl-PPMOl-PEPMl 
EPMD2-PPM02-PEPM2 
EPM03-PPM03-PEPM3 
DPMDl-PPMOl-Pl 
, OPM02-PPH02-P2 
0PMD3-PPM03-P3 

INCREMENT COUNTER FOR NUMBER OF ITERATIONS FOR CONVERGENCE 
J-2 

IF (IROBUST.EQ.Z) J-5 
CALL KTITER(K*J) 

C . 

C 
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HY TAYLOR-SERIES APPROXIMATED POSTERIOR MEAN AND COVARIANCE 
MATRICES 

100 PAPHl-PEPNl 
PAPM2-PEPM2 
PAPH3-PEPM3 
K-1 

lOS PLI-PAPMl 
PL2«PAPM2 
PL3-PAPH3 

PAPM12-PAPM14PAPH2 
PAPM13*PAPM1-»PAPH3 
PAPM23-PAPM24PAPM3 
IF CPAPM12-1.E-1A) 107.107.108 
lOT TEMP-0. 

60 TO 10<) 

108 TEMP-212/PAPM12 

109 PAPM1-(ZIN*PAPM1A(TEMP*Z13/PAPM131 WSSN 
PAPM2-(Z2N^PAPM2*<TEMP*Z23/PAPM23»»/SSN 
PAPM3-1.-PAPN1-PAPM2 

IF (PAPM3-0.1 121.11A.11A 

114 IF *PAP«l-0.1> 115,115.116 

Its IF (ABS(PAPM1>PL1I-0«00001I 117.120.120 

116 IF (A8S<PAPM1-PL1I/PAPH1-C0MVCRII 117.120.120 

117 IF (PAPM2«0.1I 116.118.119 

118 IF (ABS(PAPM2-PL2I>0. 000011 125.120.120 

119 IF (ABS(PAPM2>PL2)/PAPMZ-C0NVCRn 12S.120.120 

120 K-K*l 

IF <K>1000I 105.105.121 

121 PRINT 122. K.XNU1.PI.PL1.P12.C0NVCRI.IR0BUST.NTS.TPI0.PAPH3.XNU2.P 
12.PAPM1.PAPM2 

122 FORMAT!// ♦ EXCESS NUMBER OF ITERATIONS FOR INCOMPLETE-DATA PAPM 
1. 15*130* XNU1-*F9,4.* Pl-*F9.6.3X*PL1-*F11.8.3XPPL2-*F11.8/P CONV 
2(R0I-AF7.5* IR0BUST-*I2* NTS"*I3* TPI0-*F6.2* PAPM3-*F6.4* XNU2-*F 
39. 4P P2-AF9.69 PAPM1-*F11.8* PAPM2-PF11.8> 

IF (K-10301 125.125.250 

CONVERGENCE FOR T.S. APPROX. POSTERIOR MEAN INCOMPLETE DATA 

125 EAPMNl-PAPMl-PEPMl 
< EAPMN2-PAPN2-PEPM2 
, EAPHN3-PAPM3-PEPM3 
OAPMNl-PAPMl-Pl 
... . 0APMN2-PAPM2-P2 
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0AFMN3-PAPH3-P3 

INCREMENT CaUNTER FOR NUMBER OF ITERATIONS FOR CONVERGENCE 

IF (IROBUST-II 127»128»129 

127 J-3 
GO TO 130 

128 

GO TO 130 

129 J«6 

130 CALL KTITER(K»J) 

IF (IROBUST.GT-.Ot RETURN 

APPROXIMATED POSTERIOR VAR/COV MATRIX. NONITERATIVE METHOD. 

150 P12-PAPM14PAPM2 
P13-PAPM14PAPM3 
P23«PAPM2^PAPM3 

CAUTION. INSURE THAT P12> P13t AND P23 ARE NOT IN COMMON PROM 
GENERATED Pl« P2» AND P3 

R12-PAPM1/P12 
R13-PAPM17P13 
R21-PAPM2/P12 
R23-PAPM2/P23 

SSN > SUM OF DATA PLUS SUM OF PRIOR PARAMETERS XNUI 

* » 

T«SSN*ISSNtl.) 

P12S0-P12*P12 
P13S0-P13*P13 
P23SQ»P23*P23 
ZRP21»Z12*R21/P12 
ZRP13*213PR13/P13 
ZRP12*Z12*R12/P12 
ZRP23«Z23*R23/P23 
C CALCULATE A(l»l) 

A112»(ZRP21*Z13/P13)P*27T 
A113»(ZRP21*R21I/(P12*TJ 
A11A»Z13/<P13SQ*TI 
Aa»l>«-1.+A112-A113-A11A 
C CALCULATE A(l»2) 

A121«(2RP13-ZRP12)P<ZRP2l*Z13/P13l 

A122-ZRP12PR21/P12 
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A123-ZRP13/P13 

.. AU»2I-2.*IA121*A122-A123I/T 

TENP-Z12/P12SO 
C CALCULATE A(l>3) 

A131«IZRP12-ZRP13I**2 
.. . A132*ZRP12*R12/P12^ZRP13*R13/P13 
A(1«3}-(A131-A132)/T 
C CALCULATE 

Bfl»ll — (SSNAPAPM1«P23«ZRP12«PAPH2*ZRP13*PAPN3>/T 
CALCULATE AC2»1I 

A211>(ZRP21«Z13/P13I*(ZRP23-ZRP21) 
A212»TEHP*R21**2 
A(2«li-(A211>A212I/T 
C CALCULATE A(2*2} 

A2221»TEMP*PAPM1>Z23*C1.-2.*PAPM1)/P23S0 

A2222*TEMP*PAPM2^Z13*(1.-2.^PAPH2>/P13SQ 

A222«A2221«A2222 

A223«TEMP*«Z12-2.)*PAPH1*PAPM2/P12S0 
A22A«Z13*Z23*(P12-2.*PAPMI*PAPH2I/IP13S(JAP23SQ> 
A(2*2I>(>T«A2224A223«A224I/T 
C CALCULATE A(2»3I 

A23I*(ZRP12>Z23/P23I*(ZRP13>ZRP12I 
A232«TEMP*RIZ*P2 
A(2*3)«(A231^A232)/T 
C CALCULATE B(2>1) 

B<2»II-PAPM1*PAPM2*(SSN4TEMP»/T 

C 

C CALCULATE A(3»l) 

A3ll«(ZRP21-ZRP23>**2 
A312»ZRP21*R21/P12 
A313-ZRP23PR23/P23 
A<3»ll-<A3ll-A312-A313l/T 
C CALCULATE A(3>2) 

A321»(-ZRP12-Z23/P23l*(ZRP2l-ZRP23l 

A322-ZRP12PR21/P12 

A323-ZRP23/P23 

A(3*2)>2.P(A321«A322>A323I/T 
C CALCULATE A(3>3) 

A332-(ZRP12«Z23/P23)**2 
A333«ZRP12*R12/P12 
. A33A-Z23/P23S0 

A*3»3»*(-T*A332-A333-A33^I/T 
C CALCULATE B(3»l) 

BC3#1)—CSSN*PAPH2*P13+ZRP21*PAPH1»ZRP23*PAPM3)/T 
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SOLVE SrSTEN A«X>B FOR X. X IS VECrOR OF COVARIANCES C11»C12«C22 

CALL NATINV<3>3»A*1«8*1*0ETERN»ISCALE>II>IV0T»IWK1 
IF (ABS(0ETERfO-5.0E-lA) 212*212»220 

212 PRINT 214* 0ETERN»K#XNUl#XNU2#PI*P2*PAPNl*PAPf12#((A(I*J)/J*l#3)#B( 
1I*1)*1-1«3) 

214 FORNATI///* SINGULAR SYSTEM. DETERM-*F18 .14*4 NUMBER OF ITERATIONS 
1 WAS4I2*4 XNU1-4F6.1.4 XNU2-4F6.1*4;P1-*F9.6»* P2-*F9.6»* PAPHl»*F 
29.6/4 PAPM2-4F9.6*4 &X«B MIO>C ALC I S43E20.8.5X* 4X 14* 3X* 4-«* £23.8/ 
3a4X*3E20.8*8X»4X2**3X»4««*E23.8/34X*3E20.8»5X*4X34*3X>4«**E23.8/) 
60 TO 230 

DIFFERENCE OF APPROXIMATED POSTERIOR COVARIANCES FROM EXACT 
POSTERIOR COVARIANCES. Cll* C12* AND C22 

220 APNC11«B(1*1I 
APNC12-B(2»1I 
APMC22*B(3«1) 

IF (APHC11>0.1 229*229*221 

221 IF CAPMC22-0.) 229*225*222 

222 EAPMCll-APMCll-EPMCll 
EAPMC12-APMC12-EPMC12 
EAPMC22-APMC22-EPMC22 
RETURN 

229 PRINT 226* APMCll* APMC12*APMC22*EPMC11*EPMC12*EPMC22*XNU1*XNU2»PAP 
1M1*PAPM2»NTS 

226 F0RMATI//4 APPROXIMATED VARIANCE IS NEGATIVE. APMC11>4E21.14*4 AP 
1MC12«4E18.11*4 APMC22-4E18.11/4 EPMC11-4E18. 11*4 EPMC 12»*E18 . 11 *4 
2EPMC22»4E20.13*4 XNU1»*F4.0*4 XNU2«4F4.0*4 PAPM1>4F9.3* 4 PAPM2-4F5 
4.3*4 NTS-4I3/) 

230 COVSKIP*! 

RETURN 

290 1ST0P>1 
RETURN 
END 
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SUBROUTINE CQUNTS(BIAS«RELOIFF» J) 

COUNT NUMBER OF NUNXZ UNTERNINATEO TRIALS THAT HAVE NE6ETAIVE* 
ZERO* AND POSITIVE BIAS ANO THAT HAVE ABSOLUTE RELATIVE 
DIFFERENCES LESS THAN CERTAIN PERCENTAGES. 

BIAS • (APPROX>EXACT) OR (APPROX-GENERATED P) 

RELOIFF • ABS(BIAS/EXACTt OR ABStBIAS/GENERATEO P) . 

(RECALL THAT CQV IS NEC SO WANT DENOMINATOR INCLUDED IN ABS VALUE) 

J DENOTES* IN SUBSEQUENT ORDER* ONE OF EAPMCll* EAPHC12* EAPMC22* 
(BIAS OF APPROX T.S. EXPANSION FOR EXACT POSTERIOR COV) 

DHLCl ANO 0MLC2 (COHPLETE-OATA MLE BIAS FROM GENERATED OR GIVEN P) 
(THUS* J>3 REFERS TO EAPNC22) 

COMMON DEP(3*3)*OQL(4*3)*E(A*3)*IROBUST*NTS*Pl*P2*P3*TPID 
COMMON/BIASRO/COUNTB(3*8)*COUNTRO(8*8) 

A8«ASS(8IAS) 

IF (AB-l.E-15) 3*3*1 

1 IF (BIAS-0.) 2*2*% 

C NEGATIVE BIAS 

2 COUNTB(l*J)-COUNTB(l*J)4l. 

GO TO 9 

C ZERO BIAS (CDC 6600 COMPUTER ACCURACY IS 1% SIGN FIGURES BUT 

C CONSIDER ONLY 15 DECIMAL PLACES FOR ZERO BIAS 

3 C0UNTB(2* JI«C0UNTB(2»J)«I. 

GO TO 5 

C POSITIVE BIAS 

% C0UNTB(3* J)«C0UNTB(3*J)«1. 

C 

C 25 PERCENT RELATIVE DIFFERENCE 

9 IF (RELDIFF-0.25) 8* 8*30 

8 C0UNTRD(8*J)-C0UNTR0(8*J)4^1. 

C 20 PERCENT RELATIVE DIFFERENCE 

IF (RELDIFF-0.20) 10*10*30 
10 C0UNTR0(7*J)-CQUNTR0(7*J)^1. 

C 15 PERCENT RELATIVE DIFFERENCE 

IF (RELDIFF-0.15) 12*12*30 
12 C0UNTR0(6* J)«C0UNTRD(6*J)«^1. 

C 10 PERCENT RELATIVE DIFFERENCE 

IF (RELDIFF-0.10) 1%*1%*30 
1% C0UNTR0(5*J)»C0UNTRD(5*J)»1. 

C 5 PERCENT' RELATIVE DIFFERENCE 

IF (RELDIFF-0.05) 16*16*30 
16 C0UNTRD(%*J)-C0UNTR0(%*J)«1. 
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C 1 PERCENT RELATIVE DIFFERENCE 

XF (RELOIFF-O.Oll 18>18«30 
18 COUNTROO* J)-COUNTRO(3*J)«1. 

C .1 PERCENT RELATIVE DIFFERENCE 

IF CRELOIFF-0.0011 20»20»30 
20 CDUNTRD(2f J)-CQUNTR0(2»JI^1. 

C .01 PERCENT RELATIVE DIFFERENCE 

IF (RELDIFF-.OOOl) 22»22»30 
22 COUNTRDd* JI-CQIJNTR0(1»JH>1. 

30 CONTINUE 

IF (REL0IFF»0.15) A0»31*31 

31 IF (J>AI 33.A0.32 

32 IF U-91 A0«A0»33 

33 PRINT 35* J*NTS»BIAS*REL0IFF»TPI0*IR0BUST*E(l>l)»E(l»2)#Ef4»ll»E(4 
1*21 

35 FORMATC* SUBR COUNTS. J«*I2* NTS«*I3* BIAS»*F9.7* REL0IFF«*F5. 

12P TPI0-AF6.2* IR0BUST«*I2* PEPHl* 2-*2F7.AA PAPN1*2>«2F7.A> 

40 CONTINUE 
RETUKN 
END 
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SUBROUTINE ESTNSEI Y* T2» XY«X»X2* rXHSE»N»NSE I 

CALCULATE ESTIMATES OF NSE AND SAMPLE VARIANCE OF THESE ESTIMATES 
REAL NSE(6I»MSECV 

FOR TERM-IPEI-P1M*2*<PE2-P2I**2»«PE3-P3»**2 AND CONTROL-VARIATE 
TERNCV«<PMLEC01-P1»**2>IPMLEC02-P2)**2^IPMLECD3-P3»**2 
FOR PE- DENOTING ONE OF ESTIMATORS EPM> APM* PM0» AND MLE AND 
PMLECO- DENOTING COMPLETE-DATA MLE 

Y ■ SUM OF N TERM 
Y2 • SUM OF N TERM^TERM 
XY ■ SUM OF N TERM*T£RMCV 
, X > SUM OF N TERMCV 

X2 • SUM OF N TERMCV*TERMCV 
N • NUMBER OF TERMS 

TXMSE • TRUE MEAN SO ERROR OF CONTROL VARIATE 
MSECV • USUAL SAMPLE MSE FOR THE CONTROL VARIATE 

MSE(l) IS USUAL MSE (B-0 IN MSE REGRESSION ESTIMATEI 
MSE(2) IS VAR OF MSE(l) 

MSEC3) IS USUAL CONTROL-VARIATE MSE fB«l IN MSE REGRESSION EST) 
MSEfAl IS VAR OF MSE(3I 

MSEI5I IS LEAST-SQUARES REGRESSION ESTIMATE MSE CB-LEAST-SQS ESTI 
MSEI6I IS VAR OF MSE<5I 

NOTE THAT MSE(S) SHOULD HAVE SMALLEST VARIANCE. HOWEVER* IT WILL 
BE A BIASED ESTIMATE. HENCE* USE IT IN ANALYSES ONLY IF IT 
DIFFERS FROM EITHER ONE OF TWO UNBIASED ESTIMATES BY NO MORE 
THAN IX. (IE* IT CAN DIFFER BY MORE THAN IX FROM EITHER MSEIl) OR 
MSECS) BUT NOT BOTH.) 

GENERAL FORM OF ESTIMATED MSE IS 

MSE«MSE(l)*B*(TXMSE-MSECV) 


USUAL MSE (B«0) 

XN-N^CN-1) 

HSE(1)-Y/N 

MSE(2)«(Y2-N*MSE(1)AMSE(1))/XN 
USUAL CONTROL-VARIATE MSE CB-1) 
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XM»N*<N-2I 

T0»2.*xr 

HSECV-X/N 

Tl«N*MSECV 

72»2.»M*TX»1SE A 

T3-T2*TXNSE/2. ^ 

D-TXMSE-MSECV 

«SEI3)-«SE(1»*0 

HSE(4)«(Y2-T0*X2«T2*(NSEai-nSECVH>T3-N*nSEf3l*NSE(3l)/XN 
LEAST-SQUARES REGRESS10N*ESTINATE NSE (B*LEAST-SQUARES ESTIMATE) 
B«tXY-Tl*HSE(ll)/<X2-Tl*MSECV) 

MSE(5)-MSe<l)>B*0 

MSE(6l«(Y2-6*T0*B2*X2*B*T2*tMSE(l)-B*MSECV)*B2*T3-M*MSEt5l»MSEt5)) 

1/XN 

RETURN 

END 
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SUBROUTINE KTITER(K»J) 

INCREMENT COUNTERS EOR ID AVERAGING NUMBER OF ITERATIONS AN 
ESTIMATOR REQUIRED AND (2) DETERMINING HOW MANY CASES IN A 
REPLICATION TOOK A SPECIFIED NUMBER OF ITERATIONS 

V 

K IS NUMBER OF ITERATIONS REQUIRED TO MEET CONVERGENCE CRITERION 
J DENOTES* IN SUBSEQUENT ORDER* ONE OF ESTIMATORS MLE* PMORO* 
APHRO* APMRl* PMOR2* APMR2 ITHUS* J*A REFERS TO APMRD 

COMMON 0EPI3»3I»0QLU«3)*E(A*3»*IR0BUST*NTS*P1*P2*P3*TPI0 

C0MM0N/0ATA/XDATA(2I*Z0ATA(6I 

COMMON/ ITKT/AVNUMIT 1 61 *CTNUMIT( 6*101 

FOR AVERAGING NUMBER OF ITERATIONS FOR JTH ESTIMATOR 

. AVNUMITUI-AVNUMITI JDK 

INCREMENT COUNTER FOR NUMBER OF ITERATIONS 

1*1 

IF (K>D 20*20*2 
2 1-2 

IF IK-21 20*20*3 

3 1-3 

IF (K-31 20*20*6 

♦ 1-6 

IF IK-6) 20*20*5 

5 1-5 

IF IK-5) 20*20*6 

6 1-6 

IF IK-6) 20*20*7 

7 1-7 

IF IK-7) 20*20*8 

8 1-8 

IF IK-10) 20*20*9 
9 1-9 

IF IK-15) 20*20*10 
10 1-10 

20 CTNUMITIJ*I)-CTNUMITIJ*I)n 
IF IK-25) 30*25*25 

25 PRINT 27* NTS* K* J* TPID* XDATAI1)*X0ATA 12 )* IROBUST* I IE I II * J J) » J J-l*2 
2)*II-1*6)*{ZDATAIII)*II-1*6) 

27 FORMAT!* SUBR KTITER. NTS-*I3* NUMBER OF ITER IS*I6* FOR METHOD 
1 4-*Il* IMLE*PMORO*APMRO*APMR1*PMOR2*APMR2>. TPI0-*F6.2* XIC.O.) 
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2«*2F7.2/5X* IROBUST«*Il* PEPM1*2«*2F6.4* PHLE1»2-*2F6.^« PPHOl 

32»*2F6.4* PAPH1*2>«2F6.4* 2-*6F^.0> 

SO CONTINUE 
RETURN 
END 
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SUBROUTINE SUNNARr(X»N>LE#LH*H«UH#UE) 

TUKEY"S FIVE-POINT DATA SUHHARY. ROUTINE SORTS INPUT VECTOR X OF 
LENGTH N AND THEN CALCULATES LOWER EXTREME LE* LOWER HINGE LH> 
NEOIAN N* UPPER HINGE UH> AND UPPER EXTREME UE. 

REFERENCE "EXPLORATORY DATA ANALYSIS" BY JOHN W. TUKEY 

FORTRAN EXTENDED VERSION 4.6> CDC 6600 COMPUTER (16 SIGN FIG S.P.) 

PROGRAMER IS KAREN R CREOEUR* NASA* LANGLEY RESEARCH CENTER 

DIMENSION X(NI 
REAL LE*LH*N 

SORT DATA IN ASCENDING ORDER 

CALL AS0RT(X»1«NI 

NEOIAN 

XN«(N4l*)/2. 

, . K»XN*l.E-12 

M"XIK> 

IF (ABS(XN-KI-l.E-8) 5*5*1 
1 N"(NtX(K*in/2. 

HINGES 

5 XN«(K*l)/2. 

K»XN*1.E-12 
. LH«X(KI 

UH-X(N*l-K» 

IF (ABS(XN-KI-l.E-8» 15*15*10 
. 10 LH»(LH*-X(K*-in/2. 

UH«(X(N-K)*UH)/2. 

EXTREMES 

15 LE-X(1I 
UE-X(NI 
RETURN 
END 
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SUBROUTINE BESTESTCBESTOL) 

BY TWO DIFFERENT CRITERION (SUNNED ABSOLUTE RELATIVE DIFFERENCE 
RELD— AND SUNNED SQUARED ERROR SE — FOR SUN BEING OVER THE THREE 
CQNPONENTS OF AN ESTINATOR) DETERNINE WHICH ESTINATOR IS BEST FOR 
A GIVEN ONE OF THE TRINONIAL SINULATION TRIALS 

TIES* SCORE AS BEST EACH ESTINATOR THAT TIES FOR BEST. 


THE FOUR ESTINATORS IN E SHOULD BE IN THE FOLLOWING ORDER PERN* 
PNLE* PPND* PAPN. BIASES SHOULD BE IN CORRESPONDING ORDER. 

DINENSION 6ESTOL(4»2)»ROEP(3)*RDQL(A).RELDEP(3*B)»RELDQL(A*3) 
DINENSION SEEP(3l»SEOL(AI»V(A)»W(A>»X(3l»Y(3) 

CONNON DEP(3»3)>00L(4>3l>E(A>3)»IR0aUST»NTS>Pl«P2>P3*TPI0 
CaNN0N/BEST/BESTEP(3.2)>CTR0EP*CTRDQL(3).PRDEP(9.3)>PRDQL(9.7)»SBI 
1ASEP(3«3I»SBIAS0L(3»7I 
C0NN0N/DATA/XDATA(2I»Z0ATA(6) 

IR«IR08UST>1 

IF (IROaUST-ll 1*100*150 
FOR EPN CONPARISONS 

1 DO 10 1-1»3 
$EEPm-0. 

DO 2 J<il*3 

X(I)«SEEP(I l-SEEP<II>OEP(I* J)*DEP(I*J). 

2 CONTINUE 

INCORPORATE SUBROUTINE COUNTS TWICE (ONCE FOR EACH OF EP AND QL 
CONPARISONS) IN THIS SUBROUTINE TO SAVE PROGRAN EXECUTION COST OF 
NANY SUBROUTINE CALLS AND INDEX RESETTINGS. 

DETERNINE SIGN OF FIRST CONPONENT OF ESTINATOR 

IF (ABS(DEP(I*in-l.E-13) 5*5*3 

3 IF (OEP(I»1)-0.) 4*4*5 

C NEGATIVE BIAS 

4 SBIASEP(1*I)-SBIASEP(1*I)41. 

GO TO 10 

C ZERO BIAS (CONSIDER ONLY 13 OECINAL PLACES) 

5 SBIASEP(2*I)"SBIASEP(2*I)'»1. 

GO TO 10 
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C POSITIVE BIAS 

6 SBIASEP(3>I)-S81ASEPC3»im. 

10 CONTINUE 

CALL ASQRT(X<1*3) 

DO 20 1-1*3 

IF (SEEP(n.EQ.X(l)l BESTEPf I*1)-BESTEPCI»1)41. 
kO CONTINUE 

OETERHINE WHETHER ANY PEPN COMPONENT IS ZERO 
ICK-0 

IF (E(l*ll-1.E-10) 30*25*25 

25 IF (E(1«2I-1.E«10I 30*26*26 

26 IF (E(1*3I-1.E-10I 30*32*32 
30 ICK-l 

32 IF (lCK-0» 33*33*56 

33 CTR0EP-CTRDEP>1, 

DO 35 1-1*3 
R0£Pm-0. 

DO 3A J-l*3 

RELOEP(I»J)-ABS(OEPn*JH/Efl»J) 
y<H-ROEPm-ROEP(I>*R£LOEPII>J) 

36 CONTINUE 
35 CONTINUE 

CALL AS0RT<Y*1*3> 

DO 60 1-1*3 

IF (ROEPdl.EO.Yd)) BESTEPd*2)-BESTEPd*2)«l. 

60 CONTINUE 

FOR OETERNINING PRaPORTION OF CASES FOR WHICH X ABSOLUTE RELATIVE 
DIFFERENCE FOR EACH OF ALL THREE ESTIMATOR COMPONENTS IS LESS THAN 
SPECIFIED AMOUNTS (INCORPORATED IN PART FROM SUBROUTINE COUNTS} 

DO 55 1-1*3 
II-l 

DO 53 J-l*3 

60 TO (61*63*65*67*69*51*520*528*53) II 

61 IF (RELOEPd* J)-0.0001) 53*62*62 

62 II-2 

63 IF (RELOEPd* Jl-O.OOl) 53*66*66 
66 II-3 

65 IF (REL0EP(I*J)-0.01) 53*66*66 

66 II-6 

67 IF (REL0EP(I*JI-0.05) 53*68*68 

68 II-5 
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49 IF (RELOEFd* JI-O.IOI 53t50*50 

50 Il*6 

51 IF (RELOEP(I*J)-O.I5I 53*S2»52 

52 11-7 

S20 IF (RELOEP(I>JI-0.20I 53*525»S25 
S25 II-8 

S28 IF (RELOEP(I»J}-O.23l S3»530>530 
S30 II-9 

53 CONTINUE 

PROEP(II*n«PROEP(II*im. 

SS CONTINUE 

IF fREL0EPI3>l)-0.15) S41*543>543 
541 IF (REL0EP(3>2)-0.15) 56*543»543 

543 PRINT 544* NTS* IROBUST* TPID*RELOEP f 3*1)» RELDEP (3*2 ) * ( I E f I* J ) » J-l* 2 
1)*I-1*4*3I* ( (E(I*JI*J-l*2l*I-2*3)*(f0EP(I*Ji»J>l*2)*I«l*3)*X0ATA(l 
1I*X0ATA(21* (Z0ATA(II*I«1*6) 

544 FOR«AT(4 SUBR BESTEST. NTS-4I3* IROBUST»*Il* TPI0-*F4.2* RELOE 

1PC3*1-2I«*2F5.2* PEPM1*2«42F6,4* PAPN1»2-*2F6.4* PMLE1*2»*2F5.3 

2/15X*PPN01*2«*2F7.4* OEP«P3C2F7,4»3X)** X»Z»*2F4.0*2X»6F4.0) 

FOR OL CONPARISONS* ROUTINE IS CALLED FOR EACH OF THREE 
ROBUSTNESS SETS 


ROBUSTNESS SET 0. ORIGINAL PRIOR. 

56 Il«l 
L-2 
ICK«0 

OETERNINE WHETHER ANT P CONPONENT IS ZERO 

IF <P1>1.E>10) 59,57*57 

57 IF fP2-l.E-10) 59,58,58 

58 IF CP3-1.E-10I 59*61*61 

59 ICK-1 

.61 00 75 1-11*4 
SEQL(I)-0. 

DO 610 J-1,3 

Vm-SEOLU l-SEQL(II«OQLfI* J)*DOL(I*J) 

610 CONTINUE 

IF (IR-2) 62*615*64 
615 IF (I-3I 75,75,63 
62 K-I 

GO TO 65 
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63 K-5 

60 TO 65 

66 K«I*3 

65 IF (ABS(0QL(I*lt)*l.E-13l 69*69*67 

67 IF (OQL(I*n-0.} 68*66*70 
C NEGATIVE BIAS 

68 SBIASQL(l*Kl-SBIASQLCl*Km. 

60 TO 75 

C ZERO BIAS (CONSIDER ONLY 13 OECINAL PLACES) 

69 SBtASQL(2*KI-S8lASQL(2»K)4-I. 

60 TO 75 

C POSITIVE BIAS 

70 SBIASQL(3*K)-SBIASQL(3*Kin* 

75 CONTINUE 
CALL AS0RT(V*L*6) 

DO 76 I-L*6 

IF (SEOL(I).EQ.V(L)) BESTQL(I*1)*BEST0L(I*1)«1. 

76 CONTINUE 

UMSORT V FOR ROBUSTNESS SETS 

00 770 I>L*6 
V(I)-SEQL(I) 

770 CONTINUE 

IF ANY ESTINATOR IS ZERO* SKIP OIV FOR RELATIVE DIFF FDR ALL EST 

IF (ICK.GT.O) RETURN 
CTR0QL(IR)-CTRD0L(IR)^1. 

DO 78 I>I1*6 

REL0QLM*1)«ABS(DQL(I*1))/P1 
RELDQL(I*2)-ABS(DQL(I*2) )TP2 
RELOOL(I*3)>A8S(OOL(I*3))/P3 
. ROOLCD-O. 

00 77 J-l*3 

W(I)«ROQL(I)>ROQL(I)^RELOQL(I*J) 

77 CONTINUE 

78 CONTINUE 
CALL AS0RT(W*L*6) 

DO 79 I-L*6 

IF (ROQL(I).EQ.W(LI) 8ESTQL(I*2)-BESTOL(I*2)'M. 

79 CONTINUE 
C 

C UNSORT W FOR ROBUSTNESS SETS 
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OF POOR QUAlOT 

00 790 I-L.V 
• Udl-ROQLdl 
790 CONTINUE 

IF (IROSUST. E0.lt 11-4 
00 9d I-I1»4 
. ZI-1 

IF d(l*2t 80*81*82 

80 K-1 

60 TO 63 

81 K-9 

60 TO 83 

82 K-I«3 

83 00 96 J-l*3 

60 TO (84*86*86*90*92*94*950*958*96) II 

84 IF (REL0OLd*J)*0.0O01) 96*85*85 

85 II-2 

86 IF (REl00Ld*J)>0.001) 96*87*87 

87 II-3 

88 IF (REIOQUI* J)-0.Q1) 96*89*69 

89 II-4 

90 IF (REL00Ld*J)-0.05) 96*91*91 

91 11-5 

92 IF (RELOOU1*J)>0.10) 96*93*93 

93 II-6 

94 IF IRElOOLd*J)-0.15) 96*95*95 

95 II-7 V, 

950 IF (REL0QLd*JI-0.20) 96*955*955 

955 II-8 

958 IF (RELOOU 1*J)>0.25) 96*960*960 

960 II-9 

96 CONTINUE 

PR0QLdI*KI-PR00UII*K)4-l. 

98 CONTINUE 
RETURN 

ROBUSTNESS SET 1. UNIFORM PRIOR. 

100 11-3 
L-2 

SET POSTERIOR NODE EQUAL TO N.L.E. 

00 102 1-1*3 
0QL(3*I)-0QL(2*I) 

102 CONTINUE 



60 TO 61 
C 

C ROBUSTNESS SET PERTURBED PRIOR 
C 

ISO 11«3 

60 TO 61 
END 
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LANGLEY LIBRARY FUNCTION URAN 

Language : COMPASS 

Purpose : URAN generates uniformly distributed random numbers over the 

interval ( 0 , 1 ). 

Use : Y = URAN(X) 


X An input real number on which three conditions exist: 

X = 0, The next random number is generated and returned. If no 
previous call was made, a default seed of 
17171274321477413155B is provided. 

X < 0, A random number is not generated but the last previously 
generated random number (or the seed if no random number 
has been generated) is returned. 

X > 0, The exponent part of X is set to 1717B and the low order 
bit is set to one. This result is returned as the seed 
of a new sequence, and any additional calls to URAN will 
be based on a sequence using this seed. 

Method : This pseudorandom-number generator is multiplicative with 

algorithm 

X^.^j = 43490275647445 X- mod(2^®). 

Each random number is generated from the previous one by taking the 

lower order 48 bits of the 96 bit product produced by 

X ^^2 = 43490275647445 X^ . The exponent of the product is such that 

X^^^ is constrained to lie between 0 and 1 . 

46 

Accuracy : The generator has a full period of 2 . Extensive statisti- 

cal testing for randomness and distribution were performed to establish 
its validity as a reliable random number generator. 

SPECTRAL NUMBERS: C^ C 3 C^ C^ 


2.839 2.095 


1.819 0.978 
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References : (a) Ahrens, J. H. and Dieter, U.: "Computer Methods for 

Sampling from Gamma, Beta, Poisson, and Binomial Distributions," 
Computing 12, 1974, p 224. 

(b) Ahrens, J. H., Dieter, U., and Grube: "Pseudo-Random 

Numbers: A New Proposal for the Choice of Multi pi icators," Computing 

6, 1970, pp 121-138. 

(c) Knuth, Donald E.: The Art of Computer Programming, 

Vol . 2 (Semi numerical Algorithms). Addison-Wesley , Reading, Mass. 1969. 

Storage : 13 octal locations 


Subroutine date: March 1, 1977 
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LANGLEY LIBRARY SUBROUTINE URANV 


Language : COMPASS 

Purpose : URANV generates uniformly distributed random numbers over the 

interval (0,1). 

Use: CALL URANV(X,N,V) 


X An input real number on which three conditions exist: 

X = 0, A vector of random numbers is generated using the last 
random number generated on the previous call as a seed. 

If no previous call was made, a default seed of 
17171274321477413155B is provided. 

X < 0, The last random number calculated by the routine, or the 
default seed if no previous call was made, is returned in 
V(l). V(2), ..., V(N) are not altered. 

X > 0, The first random number is found by packing an exponent 
of 1717B and the coefficient' part of X into V(l), and 
setting the low order bit to one. Random numbers V(2), 
..., V(N) are then calculated using the algorithm given 
under METHOD. 

M Input integer specifying the number of random numbers to be re- 
turned in V. 

N j< 1, V(l) is calculated and returned. 

N > 1, V(l), ..., V(N) are calculated and returned. 

V An output one-dimensional real array dimensioned at least N. On 
output, V will contain the N calculated random numbers. 


Method : This pseudorandom-number generator is multiplicative with 

algori thm 

X^^j = 43490275647445 X^. mod(2^®). 

Each random number is generated from the previous one by taking the 
lower order 48 bits of the 96 bit product produced by 

= 43490275647445 X^ . The exponent of the product is such that 

X^^j^ is constrained to lie between 0 and 1. 
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46 

Accuracy : The generator has a full period of 2 . Extensive statisti- 

cal testing for randomness and distribution were performed to establish 
its validity as a reliable random number generator. 


SPECTRAL NUMBERS: 

^2 

<:3 

<=4 

S 

t \ ^ . 

2.839 

2.095 

1.819 

0.978 


References : (a) Ahrens, J. H. and Dieter, U.: "Computer Methods for 

Sampling from Gamma, Beta, Poisson, and Binomial Distributions," 
Computing 12, 1974, p 224. 


(b) Ahrens, J. H., Dieter, U., and Grube: "Pseudo-Random 
Numbers: A New Proposal for the Choice of Multi plica tors," Computing 
6, 1970, pp 121-138. 

(c) Knuth, Donald E.: The Art of Computer Programming, 
Vol . 2 (Seminumerical Algorithms). Addison-Wesley , Reading, Mass., 
1969. 


Storage : 25 octal locations 


Subroutine date: March 1, 1977 
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LANGLEY LIBRARY SUBROUTINE MATINV 


Language ; FORTRAN 

Purpose ; MATINV solves the matrix equation AX = B, where A is a square coef- 
ficient matrix and B is a matrix of constant vectors. The solution to a 
set of simultaneous equations , the matrix inverse , and the determinant may be 
obtained. 

Use : CALL MATINV(MAX,N, A, M, B, lOP, DETERM, ISCALE, IPIVOT,IWK) 


MAX 

N 


The maximum order of A as stated in the DIMENSION statement 
of the calling program 

The order of A; 1 N 1 MAX 


A 


A two-dimensional array of coefficients. On return to the 
calling program, A~^ is stored in A. 


M The number of column vectors in B. On return to the calling 

program, X is stored in B if M > 0; for M = 0, the sub- 
routine is used only for inversion. 

B A two-dimensional array of the constant vectors B. On return to 

the calling program, X is stored in B. 


lOP 

DETERM 

ISCALE 

IPIVOT 

IWK 


Option to compute the determinant: 

0 Compute the determinant. 

1 Do not compute the determinant . 

Gives the value of the determinant by the formula 

Det(a) = 10C00 ><ISCALE)(deteRM) when lOP = 0. For lOP = 1, 
the determinant is set to 1 . For a singular matrix and 
lOP =0 or lOP = 1, the determinant is set to zero. 

A scale factor computed by the subroutine to keep the results of 
computation within the floating-point word size of the computer 

A one-dimensional array of temporary storage used by the 
subroutine 

A two-dimensional array of temporary storage used by the 
subroutine 


Restrictions ; Arrays A, B, IPIVOT, and INDEX have variable dimensions in the 
subroutine. The maximum size of these arrays must be specified in a DIMENSION 
statement of the calling program as A(MAX,MAX), B(MAX,M), IPIVOT(MAX) , and 
IWK(MAX,2). The original matrices A and B are destroyed. They must be saved 
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by the user if there is further need for them. The determinant is set to zero 
for a singular matrix. 

Method ; Jordan's method is used to reduce a matrix A to the identity matrix I 
through a succession of elementary transformations S'n-1> • • • i • If 

these transformations are simultaneously applied to I and to a matrix B of con- 
stant vectors, the results are A“1 and X where AX = B. Each transformation 
is selected so that the largest element is used in the pivotal position. (See 
ref. (a).) 

Accuracy ; Total pivotal strategy is used to minimize the rounding errors; how- 
ever, the accuracy of the final results depends upon how well-conditioned the 
original matrix is. A return with DETERM i 0 does not guarantee accuracy 
in the solutions of inverse. 

Reference ; (a) Fox, L.: An Introduction to Numerical Linear Algebra. Oxford 

Univ. Press, 1965. 

Storage ; 516 octal locations 


Subroutine date; January 1, 1975 
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